Surface Fluxes’ Diurnal Cycle

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 [bucket] = bucket_hydrology (physcon, forcvar, fluxvar, bucket, dt)

% Bucket model hydrology

% ------------------------------------------------------
% Input
%   dt                    ! Time step (s)
%   physcon.mmh2o         ! Molecular mass of water (kg/mol)
%   forcvar.rain          ! Rainfall (kg H2O/m2/s)
%   forcvar.snow          ! Snowfall (kg H2O/m2/s)
%   fluxvar.etflx         ! Evapotranspiration (mol H2O/m2/s)
%   bucket.snow_melt      ! Snow melt (kg H2O/m2/s)
%   bucket.soil_water_max ! Maximum soil water (kg H2O/m2)
%
% Input/output
%   bucket.snow_water     ! Snow water (kg H2O/m2)
%   bucket.soil_water     ! Soil water (kg H2O/m2)
%
% Output
%   bucket.runoff         ! Runoff (kg H2O/m2/s)
% ------------------------------------------------------

% --- Save current water for conservation check

snow0 = bucket.snow_water;
soil0 = bucket.soil_water;
tot0 = snow0 + soil0;

% --- Update snow and soil water for snow melt

bucket.snow_water = bucket.snow_water - bucket.snow_melt * dt;
bucket.soil_water = bucket.soil_water + bucket.snow_melt * dt;

% --- Update snow and soil water for precipitation

bucket.soil_water = bucket.soil_water + forcvar.rain * dt;
bucket.snow_water = bucket.snow_water + forcvar.snow * dt;

% --- Update snow and soil water for evaporative flux

% Evaporative loss (kg H2O/m2/s)

evap = fluxvar.etflx * physcon.mmh2o;

% Apply evaporative flux to snow, but if more snow sublimates than
% is present take the extra water from the soil as evaporation

subl = min (evap, bucket.snow_water/dt);
evap = evap - subl;

% Update snow and soil

bucket.snow_water = bucket.snow_water - subl * dt;
bucket.soil_water = bucket.soil_water - evap * dt;

% --- Runoff is soil water in excess of a maximum capacity

bucket.runoff = max ((bucket.soil_water - bucket.soil_water_max)/dt, 0);
bucket.soil_water = bucket.soil_water - bucket.runoff * dt;

% --- Conservation check

% Snow water

delta = bucket.snow_water - snow0;
err = (forcvar.snow - bucket.snow_melt - subl) * dt - delta;
if (abs(err) > 1e-03)
   error ('Bucket model snow conservation error')
end

% Soil water

delta = bucket.soil_water - soil0;
err = (forcvar.rain + bucket.snow_melt - evap - bucket.runoff) * dt - delta;
if (abs(err) > 1e-03)
   error ('Bucket model soil conservation error')
end

% Total water

delta = (bucket.snow_water + bucket.soil_water) - tot0;
err = (forcvar.rain + forcvar.snow - fluxvar.etflx * physcon.mmh2o - bucket.runoff) * dt - delta;
if (abs(err) > 1e-03)
   error ('Bucket model total conservation error')
end

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] = most (physcon, forcvar, surfvar, fluxvar, x)

% Use Monin-Obukhov similarity 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.eref       ! Vapor pressure at reference height (Pa)
%   forcvar.pref       ! Atmospheric pressure (Pa)
%   forcvar.mmair      ! Molecular mass of air at reference height (kg/mol)
%   fluxvar.tsrf       ! Surface temperature (K)
%   fluxvar.esrf       ! Surface vapor pressure (Pa)
%   fluxvar.z0m        ! Roughness length for momentum (m)
%   fluxvar.z0c        ! Roughness length for scalars (m)
%   fluxvar.disp       ! Displacement height (m)
% Output
%   fluxvar.ustar      ! Friction velocity (m/s)
%   fluxvar.tstar      ! Temperature scale (K)
%   fluxvar.qstar      ! Water vapor scale (mol/mol)
%   fluxvar.obu        ! Obukhov length (m)
%   fx                 ! Change in Obukhov length (x - obu)
% -------------------------------------------------------------------------

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

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

% --- Calculate z-d at the reference height, because this is used many times

z_minus_d = forcvar.zref - fluxvar.disp;

% --- Evaluate psi for momentum at the reference height (zref-disp) and surface (z0m)

[psi_m_zref] = psi_m_monin_obukhov (z_minus_d / x);
[psi_m_z0m] = psi_m_monin_obukhov (fluxvar.z0m / x);
psim = -psi_m_zref + psi_m_z0m;

% --- Evaluate psi for scalars at the reference height (zref-disp) and surface (z0c)

[psi_c_zref] = psi_c_monin_obukhov (z_minus_d / x);
[psi_c_z0c] = psi_c_monin_obukhov (fluxvar.z0c / x);
psic = -psi_c_zref + psi_c_z0c;

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

zlog_m = log(z_minus_d / fluxvar.z0m);
zlog_c = log(z_minus_d / fluxvar.z0c);

fluxvar.ustar = forcvar.uref * physcon.vkc / (zlog_m + psim);
fluxvar.tstar = (forcvar.thref - fluxvar.tsrf) * physcon.vkc / (zlog_c + psic);
fluxvar.qstar = (forcvar.eref - fluxvar.esrf) / forcvar.pref * physcon.vkc / (zlog_c + psic);
tvstar = fluxvar.tstar + 0.61 * forcvar.thref * fluxvar.qstar * (physcon.mmh2o / forcvar.mmair);

% --- Calculate Obukhov length (m)

fluxvar.obu = fluxvar.ustar^2 * forcvar.thvref / (physcon.vkc * physcon.grav * tvstar);

% --- Calculate change in obu

fx = x - fluxvar.obu;

function [soilvar] = phase_change (physcon, soilvar, dt)

% Adjust temperatures for phase change. Freeze or melt ice using
% energy excess or deficit needed to change temperature to the
% freezing point.

% ------------------------------------------------------
% Input
%   dt                      ! Time step (s)
%   physcon.hfus            ! Heat of fusion for water at 0 C (J/kg)
%   physcon.tfrz            ! Freezing point of water (K)
%   soilvar.nsoi            ! Number of soil layers
%   soilvar.dz              ! Soil layer thickness (m)
%   soilvar.cv              ! Volumetric heat capacity (J/m3/K)
%
% Input/output
%   soilvar.tsoi            ! Soil temperature (K)
%   soilvar.h2osoi_liq      ! Unfrozen water, liquid (kg H2O/m2)
%   soilvar.h2osoi_ice      ! Frozen water, ice (kg H2O/m2)
%
% Output
%   soilvar.hfsoi           ! Soil phase change energy flux (W/m2)
% ------------------------------------------------------

% --- Initialize total soil heat of fusion to zero

soilvar.hfsoi = 0;

% --- Now loop over all soil layers to calculate phase change

for i = 1:soilvar.nsoi

   % --- Save variables prior to phase change

   wliq0 = soilvar.h2osoi_liq(i);     % Amount of liquid water before phase change
   wice0 = soilvar.h2osoi_ice(i);     % Amount of ice before phase change
   wmass0 = wliq0 + wice0;            % Amount of total water before phase change
   tsoi0 = soilvar.tsoi(i);           % Soil temperature before phase change

   % --- Identify melting or freezing layers and set temperature to freezing

   % Default condition is no phase change (imelt = 0)

   imelt = 0;

   % Melting: if ice exists above melt point, melt some to liquid.
   % Identify melting by imelt = 1

   if (soilvar.h2osoi_ice(i) > 0 & soilvar.tsoi(i) > physcon.tfrz)
      imelt = 1;
      soilvar.tsoi(i) = physcon.tfrz;
   end

   % Freezing: if liquid exists below melt point, freeze some to ice.
   % Identify freezing by imelt = 2

   if (soilvar.h2osoi_liq(i) > 0 & soilvar.tsoi(i) < physcon.tfrz)
      imelt = 2;
      soilvar.tsoi(i) = physcon.tfrz;
   end

   % --- Calculate energy for freezing or melting

   % The energy for freezing or melting (W/m2) is assessed from the energy
   % excess or deficit needed to change temperature to the freezing point.
   % This is a potential energy flux, because cannot melt more ice than is
   % present or freeze more liquid water than is present.
   %
   % heat_flux_pot > 0: freezing; heat_flux_pot < 0: melting

   if (imelt > 0)
      heat_flux_pot = (soilvar.tsoi(i) - tsoi0) * soilvar.cv(i) * soilvar.dz(i) / dt;
   else
      heat_flux_pot = 0;
   end

   % Maximum energy for melting or freezing (W/m2)

   if (imelt == 1)
      heat_flux_max = -soilvar.h2osoi_ice(i) * physcon.hfus / dt;
   end

   if (imelt == 2)
      heat_flux_max = soilvar.h2osoi_liq(i) * physcon.hfus / dt;
   end

   % --- Now freeze or melt ice

   if (imelt > 0)

      % Change in ice (kg H2O/m2/s): freeze (+) or melt (-)

      ice_flux = heat_flux_pot / physcon.hfus;

      % Update ice (kg H2O/m2)

      soilvar.h2osoi_ice(i) = wice0 + ice_flux * dt;

      % Cannot melt more ice than is present

      soilvar.h2osoi_ice(i) = max(0, soilvar.h2osoi_ice(i));

      % Ice cannot exceed total water that is present

      soilvar.h2osoi_ice(i) = min(wmass0, soilvar.h2osoi_ice(i));

      % Update liquid water (kg H2O/m2) for change in ice

      soilvar.h2osoi_liq(i) = max(0, (wmass0-soilvar.h2osoi_ice(i)));

      % Actual energy flux from phase change (W/m2). This is equal to
      % heat_flux_pot except if tried to melt too much ice.

      heat_flux = physcon.hfus * (soilvar.h2osoi_ice(i) - wice0) / dt;

      % Sum energy flux from phase change (W/m2)

      soilvar.hfsoi = soilvar.hfsoi + heat_flux;

      % Residual energy not used in phase change is added to soil temperature

      residual = heat_flux_pot - heat_flux;
      soilvar.tsoi(i) = soilvar.tsoi(i) - residual * dt / (soilvar.cv(i) * soilvar.dz(i));

      % Error check: make sure actual phase change does not exceed permissible phase change

      if (abs(heat_flux) > abs(heat_flux_max))
         error ('Soil temperature energy conservation error: phase change')
      end

      % Freezing: make sure actual phase change does not exceed permissible phase change
      % and that the change in ice does not exceed permissible change

      if (imelt == 2)

         % Energy flux (W/m2)

         constraint = min(heat_flux_pot, heat_flux_max);
         err = heat_flux - constraint;
         if (abs(err) > 1e-03)
            error ('Soil temperature energy conservation error: freezing energy flux')
         end

         % Change in ice (kg H2O/m2)

         err = (soilvar.h2osoi_ice(i) - wice0) - constraint / physcon.hfus * dt;
         if (abs(err) > 1e-03)
            error ('Soil temperature energy conservation error: freezing ice flux')
         end
      end

      % Thawing: make sure actual phase change does not exceed permissible phase change
      % and that the change in ice does not exceed permissible change

      if (imelt == 1)

         % Energy flux (W/m2)

         constraint = max(heat_flux_pot, heat_flux_max);
         err = heat_flux - constraint;
         if (abs(err) > 1e-03)
            error ('Soil temperature energy conservation error: thawing energy flux')
         end

         % Change in ice (kg H2O/m2)

         err = (soilvar.h2osoi_ice(i) - wice0) - constraint / physcon.hfus * dt;
         if (abs(err) > 1e-03)
            error ('Soil temperature energy conservation error: thawing ice flux')
         end
      end

   end

end

function [phi_c] = phi_c_monin_obukhov (x)

% --- Evaluate the Monin-Obukhov phi function for scalars at x

if (x < 0)
   phi_c = (1 - 16 * x)^(-0.5);
else
   phi_c = 1 + 5 * x;
end

function [phi_m] = phi_m_monin_obukhov (x)

% --- Evaluate the Monin-Obukhov phi function for momentum at x

if (x < 0)
   phi_m = (1 - 16 * x)^(-0.25);
else
   phi_m = 1 + 5 * x;
end

function [psi_c] = psi_c_monin_obukhov (x)

% --- Evaluate the Monin-Obukhov psi function for scalars at x

if (x < 0)
   y = (1 - 16 * x)^0.25;
   psi_c = 2 * log((1 + y^2)/2);
else
   psi_c = -5 * x;
end

function [psi_hat_c] = psi_c_rsl (z, h, L, c1, c2)

% --- Evaluate the roughness sublayer (RSL) function psi_hat for scalars
% at z. Note that z has already been adjusted for the displacement height
% (i.e., using z - d).

% ------------------------------------------------------
% Input
%   z            ! Vertical height - displacement height (m)
%   h            ! Canopy height - displacement height (m)
%   L            ! Obukhov length (m)
%   c1           ! Parameter for RSL function phi_hat (dimensionless)
%   c2           ! Parameter for RSL function phi_hat (dimensionless)
%
% Output
%   psi_hat_c    ! RSL psi_hat function for scalars (dimensionless)
% ------------------------------------------------------

% The function to integrate depends on unstable (f1) or stable (f2)

f1 = @(x) (1-16*x/L).^(-0.5) .* (1 - (1 - c1*exp(-c2*x/(2*h)))) ./ x;
f2 = @(x) (1+5*x/L)          .* (1 - (1 - c1*exp(-c2*x/(2*h)))) ./ x;

% Numerically integrate the function from z to infinity

if (L < 0)
   psi_hat_c = integral (f1, z, inf);
else
   psi_hat_c = integral (f2, z, inf);
end

function [psi_m] = psi_m_monin_obukhov (x)

% --- Evaluate the Monin-Obukhov psi function for momentum at x

if (x < 0)
   y = (1 - 16 * x)^0.25;
   psi_m = 2 * log((1 + y)/2) + log((1 + y^2)/2) - 2 * atan(y) + pi / 2;
else
   psi_m = -5 * x;
end

function [psi_hat_m] = psi_m_rsl (z, h, L, c1, c2)

% --- Evaluate the roughness sublayer (RSL) function psi_hat for momentum
% at z. Note that z has already been adjusted for the displacement height
% (i.e., using z - d).

% ------------------------------------------------------
% Input
%   z            ! Vertical height - displacement height (m)
%   h            ! Canopy height - displacement height (m)
%   L            ! Obukhov length (m)
%   c1           ! Parameter for RSL function phi_hat (dimensionless)
%   c2           ! Parameter for RSL function phi_hat (dimensionless)
%
% Output
%   psi_hat_m    ! RSL psi_hat function for momentum (dimensionless)
% ------------------------------------------------------

% The function to integrate depends on unstable (f1) or stable (f2)

f1 = @(x) (1-16*x/L).^(-0.25) .* (1 - (1 - c1*exp(-c2*x/(2*h)))) ./ x;
f2 = @(x) (1+5*x/L)           .* (1 - (1 - c1*exp(-c2*x/(2*h)))) ./ x;

% Numerically integrate the function from z to infinity

if (L < 0)
   psi_hat_m = integral (f1, z, inf);
else
   psi_hat_m = integral (f2, z, inf);
end

function [fluxvar, fx] = rsl (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.eref        ! Vapor pressure at reference height (Pa)
%   forcvar.pref        ! Atmospheric pressure (Pa)
%   forcvar.mmair       ! Molecular mass of air at reference height (kg/mol)
%   fluxvar.tsrf        ! Surface temperature (K)
%   fluxvar.esrf        ! Surface vapor pressure (Pa)
%   surfvar.Lc          ! Canopy density length scale (m)
%   surfvar.hc          ! Canopy height (m)
%   surfvar.rc          ! Leaf Nusselt number (heat) or Stanton number (scalar)
%
% Output
%   fluxvar.z0m         ! Roughness length for momentum (m)
%   fluxvar.z0c         ! Roughness length for scalars (m)
%   fluxvar.disp        ! Displacement height (m)
%   fluxvar.ustar       ! Friction velocity (m/s)
%   fluxvar.tstar       ! Temperature scale (K)
%   fluxvar.qstar       ! Water vapor scale (mol/mol)
%   fluxvar.obu         ! Obukhov length (m)
%   fx                  ! Change in Obukhov length (x - obu)
% -------------------------------------------------------------------------

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

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

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

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

beta_neutral = 0.35;

% Lc/obu

LcL = surfvar.Lc/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('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('stable case: error in beta')
   end

end

% Place limits on beta

beta_val = min(0.5, max(beta_val,0.2));

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

dp = beta_val^2 * surfvar.Lc;                  % dp = hc - disp
fluxvar.disp = max(surfvar.hc - 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 - fluxvar.disp;
h_minus_d = surfvar.hc - fluxvar.disp;

% --- Turbulent Prandlt number (Pr) at canopy height

Prn = 0.5;         % Neutral value for Pr
Prvr = 0.3;        % Magnitude of variation of Pr with stability
Prsc = 2.0;        % Scale of variation of Pr with stability

Pr = Prn + Prvr * tanh(Prsc*surfvar.Lc/x);

% --- The "f" parameter relates the length scale of the scalar (heat) to that of momentum 

f = (sqrt(1 + 4 * surfvar.rc * Pr) - 1) / 2;

% --- 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)

c2 = 0.5;

% c1 for momentum and scalars (dimensionless)

c1m = (1 -    physcon.vkc / (2 * beta_val * phi_m_hc)) * exp(c2/2);
c1c = (1 - Pr*physcon.vkc / (2 * beta_val * phi_c_hc)) * exp(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, c1m, c2);  % momentum at (zref-disp)
[psi_m_rsl_hc]   = psi_m_rsl (h_minus_d, h_minus_d, x, c1m, c2);  % momentum at (hc-disp)

[psi_c_rsl_zref] = psi_c_rsl (z_minus_d, h_minus_d, x, c1c, c2);  % scalars at (zref-disp)
[psi_c_rsl_hc]   = psi_c_rsl (h_minus_d, h_minus_d, x, c1c, 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 + physcon.vkc / beta_val * Pr / f;

fluxvar.ustar = forcvar.uref * physcon.vkc / (zlog + psim);
fluxvar.tstar = (forcvar.thref - fluxvar.tsrf) * physcon.vkc / (zlog + psic);
fluxvar.qstar = (forcvar.eref - fluxvar.esrf) / forcvar.pref * physcon.vkc / (zlog + psic);
tvstar = fluxvar.tstar + 0.61 * forcvar.thref * fluxvar.qstar * (physcon.mmh2o / forcvar.mmair);

% --- Calculate Obukhov length (m)

fluxvar.obu = fluxvar.ustar^2 * forcvar.thvref / (physcon.vkc * physcon.grav * tvstar);
fx = x - fluxvar.obu;

% --- Roughness lengths z0m and z0c (m)

% z0m - Use bisection to find z0m, which lies between aval and bval, and refine the
% estimate until the difference is less than err

aval = surfvar.hc;
bval = 0;
err = 1e-08;

[psi_m_z0m] = psi_m_monin_obukhov (aval / x);
z0m = h_minus_d * exp(-physcon.vkc/beta_val) * exp(-psi_m_hc + psi_m_z0m) * exp(psi_m_rsl_hc);
fa = z0m - aval;

[psi_m_z0m] = psi_m_monin_obukhov (bval / x);
z0m = h_minus_d * exp(-physcon.vkc/beta_val) * exp(-psi_m_hc + psi_m_z0m) * exp(psi_m_rsl_hc);
fb = z0m - bval;

if (fa * fb > 0)
   error('RSL bisection error: f(a) and f(b) do not have opposite signs')
end

while (abs(bval-aval) > err)
   cval = (aval + bval) / 2;

   [psi_m_z0m] = psi_m_monin_obukhov (cval / x);
   z0m = h_minus_d * exp(-physcon.vkc/beta_val) * exp(-psi_m_hc + psi_m_z0m) * exp(psi_m_rsl_hc);
   fc = z0m - cval;

   if (fa * fc < 0)
      bval = cval; fb = fc;
   else
      aval = cval; fa = fc;
   end
end

fluxvar.z0m = cval;

% z0c - Use bisection to find z0c, which lies between aval and bval, and refine the
% estimate until the difference is less than err

aval = surfvar.hc;
bval = 0;

[psi_c_z0c] = psi_c_monin_obukhov (aval / x);
z0c = h_minus_d * exp(-physcon.vkc/beta_val*Pr/f) * exp(-psi_c_hc + psi_c_z0c) * exp(psi_c_rsl_hc);
fa = z0c - aval;

[psi_c_z0c] = psi_c_monin_obukhov (bval / x);
z0c = h_minus_d * exp(-physcon.vkc/beta_val*Pr/f) * exp(-psi_c_hc + psi_c_z0c) * exp(psi_c_rsl_hc);
fb = z0c - bval;

if (fa * fb > 0)
   error('RSL bisection error: f(a) and f(b) do not have opposite signs')
end

while (abs(bval-aval) > err)
   cval = (aval + bval) / 2;

   [psi_c_z0c] = psi_c_monin_obukhov (cval / x);
   z0c = h_minus_d * exp(-physcon.vkc/beta_val*Pr/f) * exp(-psi_c_hc + psi_c_z0c) * exp(psi_c_rsl_hc);
   fc = z0c - cval;

   if (fa * fc < 0)
      bval = cval; fb = fc;
   else
      aval = cval; fa = fc;
   end
end

fluxvar.z0c = cval;

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 [soilvar, fluxvar, bucket] = soil_temperature (physcon, soilvar, fluxvar, bucket, dt, f0, df0)

% Use an implicit formulation with the surface boundary condition specified
% as the soil heat flux (W/m2; positive into the soil) to solve for soil
% temperatures at time n+1.
%
% Calculate soil temperatures as:
%
%      dT   d     dT 
%   cv -- = -- (k --)
%      dt   dz    dz 
%
% where: T = temperature (K)
%        t = time (s)
%        z = depth (m)
%        cv = volumetric heat capacity (J/m3/K)
%        k = thermal conductivity (W/m/K)
%
% Set up a tridiagonal system of equations to solve for dT at time n+1, 
% where the temperature equation for layer i is:
%
%   d_i = a_i*[dT_i-1] n+1 + b_i*[dT_i] n+1 + c_i*[dT_i+1] n+1
%
% The energy flux into the soil is specified as a linear function of
% the temperature of the first soil layer T(1) using the equation:
%
%   gsoi = f0([T_1] n) + df0 / dT * ([T_1] n+1 - [T_1] n)
%
% where f0 is the surface energy balance evaluated with T(1) at time n
% and df0/dT is the temperature derivative of f0.
%
% If snow is present and the surface temperature is above freezing, reset
% the temperature to freezing and use the residual energy to melt snow.

% ------------------------------------------------------
% Input
%   dt                      ! Time step (s)
%   f0                      ! Energy flux into soil (W/m2)
%   df0                     ! Temperature derivative of f0 (W/m2/K)
%   physcon.hfus            ! Heat of fusion for water at 0 C (J/kg)
%   physcon.tfrz            ! Freezing point of water (K)
%   soilvar.nsoi            ! Number of soil layers
%   soilvar.z               ! Soil depth (m)
%   soilvar.z_plus_onehalf  ! Soil depth (m) at i+1/2 interface between layers i and i+1
%   soilvar.dz              ! Soil layer thickness (m)
%   soilvar.dz_plus_onehalf ! Thickness (m) between between i and i+1
%   soilvar.tk              ! Thermal conductivity (W/m/K)
%   soilvar.cv              ! Heat capacity (J/m3/K)
%   bucket.snow_water       ! Snow water (kg H2O/m2)
%
% Input/output
%   soilvar.tsoi            ! Soil temperature (K)
%
% Output
%   fluxvar.gsno            ! Snow melt energy flux (W/m2)
%   bucket.snow_melt        ! Snow melt (kg H2O/m2/s)
% ------------------------------------------------------

% --- Save current soil temperature

for i = 1:soilvar.nsoi
   tsoi0(i) = soilvar.tsoi(i);
end

% --- Thermal conductivity at interface (W/m/K)

for i = 1:soilvar.nsoi-1
   tk_plus_onehalf(i) = soilvar.tk(i) * soilvar.tk(i+1) * (soilvar.z(i)-soilvar.z(i+1)) / ...
   (soilvar.tk(i)*(soilvar.z_plus_onehalf(i)-soilvar.z(i+1)) + soilvar.tk(i+1)*(soilvar.z(i)-soilvar.z_plus_onehalf(i)));
end

% --- Set up tridiagonal matrix

% Top soil layer

i = 1;
a(i) = 0;
c(i) = -tk_plus_onehalf(i) / soilvar.dz_plus_onehalf(i);
b(i) = soilvar.cv(i) * soilvar.dz(i) / dt - c(i) - df0;
d(i) = -tk_plus_onehalf(i) * (soilvar.tsoi(i) - soilvar.tsoi(i+1)) / soilvar.dz_plus_onehalf(i) + f0;

% Layers 2 to nsoi-1

for i = 2:soilvar.nsoi-1
   a(i) = -tk_plus_onehalf(i-1) / soilvar.dz_plus_onehalf(i-1);
   c(i) = -tk_plus_onehalf(i) / soilvar.dz_plus_onehalf(i);
   b(i) = soilvar.cv(i) * soilvar.dz(i) / dt - a(i) - c(i);
   d(i) = tk_plus_onehalf(i-1) * (soilvar.tsoi(i-1) - soilvar.tsoi(i)) / soilvar.dz_plus_onehalf(i-1) ...
        - tk_plus_onehalf(i) * (soilvar.tsoi(i) - soilvar.tsoi(i+1)) / soilvar.dz_plus_onehalf(i);
end

% Bottom soil layer

i = soilvar.nsoi;
a(i) = -tk_plus_onehalf(i-1) / soilvar.dz_plus_onehalf(i-1);
c(i) = 0;
b(i) = soilvar.cv(i) * soilvar.dz(i) / dt - a(i);
d(i) = tk_plus_onehalf(i-1) * (soilvar.tsoi(i-1) - soilvar.tsoi(i)) / soilvar.dz_plus_onehalf(i-1);

% --- Begin tridiagonal solution: forward sweep for layers N to 1

% Bottom soil layer

i = soilvar.nsoi;
e(i) = a(i) / b(i);
f(i) = d(i) / b(i);

% Layers nsoi-1 to 2

for i = soilvar.nsoi-1: -1: 2
   den = b(i) - c(i) * e(i+1);
   e(i) = a(i) / den;
   f(i) = (d(i) - c(i) * f(i+1)) / den;
end

% Complete the tridiagonal solution to get the temperature of the top soil layer

i = 1;
num = d(i) - c(i) * f(i+1);
den = b(i) - c(i) * e(i+1);
tsoi_test = soilvar.tsoi(i) + num / den;

% --- Surface temperature with adjustment for snow melt

% Potential melt rate based on temperature above freezing

potential_snow_melt = max(0, (tsoi_test - physcon.tfrz) * den / physcon.hfus);

% Maximum melt rate is the amount of snow that is present

maximum_snow_melt = bucket.snow_water / dt;

% Cannot melt more snow than is present

bucket.snow_melt = min(maximum_snow_melt, potential_snow_melt);

% Energy flux for snow melt

fluxvar.gsno = bucket.snow_melt * physcon.hfus;

% Update temperature - If there is no snow melt, tsoi(1) = tsoi_test (as above).
% While snow is melting at the potential rate, tsoi(1) = tfrz. If snow melt is
% less than the potential rate, tsoi(1) > tfrz and < tsoi_test.

soilvar.tsoi(i) = soilvar.tsoi(i) + (num - fluxvar.gsno) / den;
dtsoi(i) = soilvar.tsoi(i) - tsoi0(i);

% --- Now complete the tridiagonal solution for layers 2 to N

for i = 2:soilvar.nsoi
   dtsoi(i) = f(i) - e(i) * dtsoi(i-1);
   soilvar.tsoi(i) = soilvar.tsoi(i) + dtsoi(i);
end

function [soilvar] = soil_thermal_properties (physcon, soilvar)

% Calculate soil thermal conductivity and heat capacity

% ------------------------------------------------------
% Input
%   physcon.hfus             ! Heat of fusion for water at 0 C (J/kg)
%   physcon.tfrz             ! Freezing point of water (K)
%   physcon.tkwat            ! Thermal conductivity of water (W/m/K)
%   physcon.tkice            ! Thermal conductivity of ice (W/m/K)
%   physcon.cvwat            ! Heat capacity of water (J/m3/K)
%   physcon.cvice            ! Heat capacity of ice (J/m3/K)
%   physcon.rhowat           ! Density of water (kg/m3)
%   physcon.rhoice           ! Density of ice (kg/m3)
%   soilvar.method           ! Use excess heat or apparent heat capacity for phase change
%   soilvar.soil_texture     ! Soil texture class
%   soilvar.sand             ! Percent sand
%   soilvar.watsat           ! Volumetric soil water content at saturation (porosity)
%   soilvar.nsoi             ! Number of soil layers
%   soilvar.dz               ! Soil layer thickness (m)
%   soilvar.tsoi             ! Soil temperature (K)
%   soilvar.h2osoi_liq       ! Unfrozen water, liquid (kg H2O/m2)
%   soilvar.h2osoi_ice       ! Frozen water, ice (kg H2O/m2)
%
% Input/output
%   soilvar.tk               ! Thermal conducitivty (W/m/K)
%   soilvar.cv               ! Volumetric heat capacity (J/m3/K)
% ------------------------------------------------------

for i = 1:soilvar.nsoi

   % --- Soil texture to process

   k = soilvar.soil_texture;

   % --- Volumetric soil water and ice

   watliq = soilvar.h2osoi_liq(i) / (physcon.rhowat * soilvar.dz(i));
   watice = soilvar.h2osoi_ice(i) / (physcon.rhoice * soilvar.dz(i));

   % Fraction of total volume that is liquid water

   fliq = watliq / (watliq + watice);

   % Soil water relative to saturation

   s = min((watliq + watice)/soilvar.watsat(k), 1);

   % --- Dry thermal conductivity (W/m/K) from bulk density (kg/m3)

   bd = 2700 * (1 - soilvar.watsat(k));
   tkdry = (0.135 * bd + 64.7) / (2700 - 0.947 * bd);

   % --- Soil solids thermal conducitivty (W/m/K)

   % Thermal conductivity of quartz (W/m/K)

   tk_quartz = 7.7;

   % Quartz fraction

   quartz = soilvar.sand(k) / 100;

   % Thermal conductivity of other minerals (W/m/K)

   if (quartz > 0.2)
      tko = 2;
   else
      tko = 3;
   end

   % Thermal conductivity of soil solids (W/m/K)

   tksol = tk_quartz^quartz * tko^(1-quartz);

   % --- Saturated thermal conductivity (W/m/K) and unfrozen and frozen values

   tksat = tksol^(1-soilvar.watsat(k)) * physcon.tkwat^(fliq*soilvar.watsat(k)) * ...
      physcon.tkice^(soilvar.watsat(k)-fliq*soilvar.watsat(k));
   tksat_u = tksol^(1-soilvar.watsat(k)) * physcon.tkwat^soilvar.watsat(k);
   tksat_f = tksol^(1-soilvar.watsat(k)) * physcon.tkice^soilvar.watsat(k);

   % --- Kersten number and unfrozen and frozen values

   if (soilvar.sand(k) < 50)
      ke_u = log10(max(s,0.1)) + 1;
   else
      ke_u = 0.7 * log10(max(s,0.05)) + 1;
   end
   ke_f = s;

   if (soilvar.tsoi(i) >= physcon.tfrz)
      ke = ke_u;
   else
      ke = ke_f;
   end

   % --- Thermal conductivity (W/m/K) and unfrozen and frozen values

   soilvar.tk(i) = (tksat - tkdry) * ke + tkdry;
   tku = (tksat_u - tkdry) * ke_u + tkdry;
   tkf = (tksat_f - tkdry) * ke_f + tkdry;

   % --- Heat capacity of soil solids (J/m3/K)

   cvsol = 1.926e06;

   % --- Heat capacity (J/m3/K) and unfrozen and frozen values

   soilvar.cv(i) = (1 - soilvar.watsat(k)) * cvsol + physcon.cvwat * watliq + physcon.cvice * watice;
   cvu = (1 - soilvar.watsat(k)) * cvsol + physcon.cvwat * (watliq + watice);
   cvf = (1 - soilvar.watsat(k)) * cvsol + physcon.cvice * (watliq + watice);

   % --- Adjust heat capacity and thermal conductivity if using apparent heat capacity

   switch soilvar.method
      case 'apparent-heat-capacity'

      % Temperature range for freezing and thawing (K)

      tinc = 0.5;

      % Heat of fusion (J/m3) - This is equivalent to ql = hfus * (h2osoi_liq + h2osoi_ice) / dz

      ql = physcon.hfus * (physcon.rhowat * watliq + physcon.rhoice * watice);

      % Heat capacity and thermal conductivity

      if (soilvar.tsoi(i) > physcon.tfrz+tinc)
         soilvar.cv(i) = cvu;
         soilvar.tk(i) = tku;
      end

      if (soilvar.tsoi(i) >= physcon.tfrz-tinc & soilvar.tsoi(i) <= physcon.tfrz+tinc)
         soilvar.cv(i) = (cvf + cvu) / 2 + ql / (2 * tinc);
         soilvar.tk(i) = tkf + (tku - tkf) * (soilvar.tsoi(i) - physcon.tfrz + tinc) / (2 * tinc);
      end

      if (soilvar.tsoi(i) < physcon.tfrz-tinc)
         soilvar.cv(i) = cvf;
         soilvar.tk(i) = tkf;
      end
   end

end

function [fluxvar, soilvar, bucket] = surface_fluxes (physcon, forcvar, surfvar, soilvar, fluxvar, bucket, dt)

% Calculate soil temperatures and surface fluxes. This routine uses the
% current estimate of surface temperature and vapor pressure to solve
% for the Obukhov length. This then provides the aerodynamic conductance,
% which is used in the surface temperature and flux calculation.

% ------------------------------------------------------
% Input
%   dt                  ! Time step (s)
%   physcon.tfrz        ! Freezing point of water (K)
%   physcon.mmh2o       ! Molecular mass of water (kg/mol)
%   physcon.hvap        ! Latent heat of evaporation (J/kg)
%   physcon.hsub        ! Latent heat of sublimation (J/kg)
%   physcon.sigma       ! Stefan-Boltzmann constant (W/m2/K4)
%   forcvar.thref       ! Potential temperature at reference height (K)
%   forcvar.uref        ! Wind speed at reference height (m/s)
%   forcvar.eref        ! Vapor pressure at reference height (Pa)
%   forcvar.pref        ! Atmospheric pressure (Pa)
%   forcvar.cpair       ! Specific heat of air at constant pressure, at reference height (J/mol/K)
%   forcvar.rhomol      ! Molar density at reference height (mol/m3)
%   surfvar.emiss       ! Surface emissivity
%   surfvar.gcan        ! Canopy conductance (mol/m2/s)
%   soilvar.method      ! Use excess heat or apparent heat capacity for phase change
%   soilvar.nsoi        ! Number of soil layers
%   soilvar.dz          ! Soil layer thickness (m)
%   soilvar.cv          ! Heat capacity (J/m3/K)
%   fluxvar.profiles    ! Use MOST or RSL for flux-profiles
%   fluxvar.qa          ! Radiative forcing (W/m2)
%   fluxvar.bucket      ! Use bucket model hydrology soil wetness factor
%   bucket.snow_water   ! Snow water (kg H2O/m2)
%   bucket.soil_water   ! Soil water (kg H2O/m2)
%   bucket.soil_water_max ! Maximum soil water (kg H2O/m2)
%   bucket.soil_beta_max  ! Soil water at which soil_beta = 1 (fraction of soil_water_max)
%
% Input/output
%   fluxvar.tsrf        ! Surface temperature (K)
%   fluxvar.esrf        ! Surface vapor pressure (Pa)
%   soilvar.tsoi        ! Soil temperature (K)
%   soilvar.h2osoi_liq  ! Unfrozen water, liquid (kg H2O/m2)
%   soilvar.h2osoi_ice  ! Frozen water, ice (kg H2O/m2)
%
% Output
%   fluxvar.rnet        ! Net radiation (W/m2)
%   fluxvar.lwrad       ! Emitted longwave radiation (W/m2)
%   fluxvar.shflx       ! Sensible heat flux (W/m2)
%   fluxvar.lhflx       ! Latent heat flux (W/m2)
%   fluxvar.etflx       ! Evapotranspiration (mol H2O/m2/s)
%   fluxvar.gsoi        ! Soil energy flux (W/m2)
%   fluxvar.gsno        ! Snow melt energy flux (W/m2)
%   fluxvar.gam         ! Aerodynamic conductance for momentum (mol/m2/s)
%   fluxvar.gac         ! Aerodynamic conductance for scalars (mol/m2/s)
%   soilvar.hfsoi       ! Soil phase change energy flux (W/m2)
%   bucket.soil_beta    ! Soil wetness factor for evapotranspiration (-)
%   bucket.snow_melt    ! Snow melt (kg H2O/m2/s)
%
% Output from Obukhov length calculation
%   fluxvar.ustar       ! Friction velocity (m/s)
%   fluxvar.tstar       ! Temperature scale (K)
%   fluxvar.qstar       ! Water vapor scale (mol/mol)
%   fluxvar.obu         ! Obukhov length (m)
%   fluxvar.z0m         ! RSL only: Roughness length for momentum (m)
%   fluxvar.z0c         ! RSL only: Roughness length for scalars (m)
%   fluxvar.disp        ! RSL only: Displacement height (m)
% ------------------------------------------------------

% --- Calculate the Obukhov length

% Calculate the Obukhov length (obu) for the current surface temperature
% and surface vapor pressure using Monin-Obukhov similarity theory or
% Harman & Finnigan (2007, 2008) roughness sublayer (RSL) theory. Use the
% functions "most" or "rsl" 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)

switch fluxvar.profiles
   case 'MOST'                     % Use Monin-Obukhov similarity theory
   func_name = 'most';             % The function name is "most", in the file most.m
   case 'RSL'                      % Use canopy coupling with roughness sublayer theory
   func_name = 'rsl';              % The function name is "rsl", in the file rsl.m
end

% Solve for the Obukhov length

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

% Uncomment this line to use MOST or RSL for neutral conditions
% [fluxvar, dummy] = most (physcon, forcvar, surfvar, fluxvar, -inf);
% [fluxvar, dummy] = rsl (physcon, forcvar, surfvar, fluxvar, -inf);

% --- Aerodynamic conductances for momentum (gam) and scalars (gac) (mol/m2/s)

fluxvar.gam = forcvar.rhomol * fluxvar.ustar * fluxvar.ustar / forcvar.uref;
fluxvar.gac = forcvar.rhomol * fluxvar.ustar * fluxvar.tstar / (forcvar.thref - fluxvar.tsrf);

% --- Calculate soil temperatures

% Save current temperatures

for i = 1:soilvar.nsoi
   tsoi0(i) = soilvar.tsoi(i);
end

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

[esat, desat] = satvap (fluxvar.tsrf-physcon.tfrz);

% Latent heat of vaporization or sublimation (J/mol)

if (bucket.snow_water > 0)
   lambda_val = physcon.hsub;              % Latent heat of sublimation (J/kg)
else
   lambda_val = physcon.hvap;              % Latent heat of vaporization (J/kg)
end
lambda_val = lambda_val * physcon.mmh2o;   % J/kg -> J/mol

% Surface conductance for water vapor (mol/m2/s)

gw = 1 / (1/surfvar.gcan + 1/fluxvar.gac);

% Soil wetness factor for evapotranspiration

switch fluxvar.bucket
   case 'no_bucket'                % No bucket hydrology
   bucket.soil_beta = 1;
   case 'use_bucket'               % Use bucket hydrology
   bucket.soil_beta = min(bucket.soil_water/(bucket.soil_beta_max*bucket.soil_water_max), 1);
end

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

fluxvar.lwrad = surfvar.emiss * physcon.sigma * fluxvar.tsrf^4;
dlwrad = 4 * surfvar.emiss * physcon.sigma * fluxvar.tsrf^3;

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

fluxvar.shflx = forcvar.cpair * (fluxvar.tsrf - forcvar.thref) * fluxvar.gac;
dshflx = forcvar.cpair * fluxvar.gac;

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

fluxvar.lhflx = lambda_val / forcvar.pref * (esat - forcvar.eref) * gw * bucket.soil_beta;
dlhflx = lambda_val / forcvar.pref * desat * gw * bucket.soil_beta;

% Net energy flux into soil (W/m2) and temperature derivative (W/m2/K)

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

% Update soil temperatures

[soilvar, fluxvar, bucket] = soil_temperature (physcon, soilvar, fluxvar, bucket, dt, f0, df0);

% --- Update surface fluxes for the change in surface temperature

dtsrf = soilvar.tsoi(1) - tsoi0(1);
fluxvar.lwrad = fluxvar.lwrad + dlwrad * dtsrf;
fluxvar.shflx = fluxvar.shflx + dshflx * dtsrf;
fluxvar.lhflx = fluxvar.lhflx + dlhflx * dtsrf;
fluxvar.gsoi = f0 + df0 * dtsrf;

% --- Adjust heat flux into soil for snow melt

fluxvar.gsoi = fluxvar.gsoi - fluxvar.gsno;

% --- Net radiation

fluxvar.rnet = fluxvar.qa - fluxvar.lwrad;

% --- Error check

err = fluxvar.rnet - fluxvar.shflx - fluxvar.lhflx - fluxvar.gsoi - fluxvar.gsno;
if (abs(err) > 1e-06)
   fprintf('err = %15.3f\n',err)
   fprintf('qa = %15.3f\n',fluxvar.qa)
   fprintf('lwrad = %15.3f\n',fluxvar.lwrad)
   fprintf('sh = %15.3f\n',fluxvar.shflx)
   fprintf('lh = %15.3f\n',fluxvar.lhflx)
   fprintf('gsoi = %15.3f\n',fluxvar.gsoi)
   fprintf('gsno = %15.3f\n',fluxvar.gsno)
   error ('surface temperature error')
end

% --- Evapotranspiration (mol H2O/m2/s)

fluxvar.etflx = fluxvar.lhflx / lambda_val;

% --- Surface vapor pressure is diagnosed from evaporative flux

fluxvar.esrf = (forcvar.eref / forcvar.pref) + fluxvar.etflx / fluxvar.gac; % mol/mol
fluxvar.esrf = fluxvar.esrf * forcvar.pref;                                 % mol/mol -> Pa

% --- Phase change for soil layers undergoing freezing of thawing

switch soilvar.method

   case 'apparent-heat-capacity'

   % No explicit phase change energy flux. This is included in the heat capacity.

   soilvar.hfsoi = 0;

   case 'excess-heat'

   % Adjust temperatures for phase change. Freeze or melt ice using energy
   % excess or deficit needed to change temperature to the freezing point.
   % The variable hfsoi is returned as the energy flux from phase change (W/m2).

   [soilvar] = phase_change (physcon, soilvar, dt);

end

% --- Check for energy conservation

% Sum change in energy (W/m2)

edif = 0;
for i = 1:soilvar.nsoi
   edif = edif + soilvar.cv(i) * soilvar.dz(i) * (soilvar.tsoi(i) - tsoi0(i)) / dt;
end

% Error check

err = edif - fluxvar.gsoi - soilvar.hfsoi;
if (abs(err) > 1e-03)
   error ('Soil temperature energy conservation error')
end

% --- Surface temperature is the first soil layer

fluxvar.tsrf = soilvar.tsoi(1);