sunML.cpp

Go to the documentation of this file.

#include "sunML.h"
 
double thin_layer_absorption(
    double R,   // dimensionless
    double T,   // dimensionless
    double I_0  // Light units such as `micromol / m^2 / s` or `J / m^2 / s`
)
{
    return I_0 * (1 - R - T);  // same units as `I_0`
}
 
double thick_layer_absorption(
    double R,   // dimensionless
    double T,   // dimensionless
    double I_0  // Light units such as `micromol / m^2 / s` or `J / m^2 / s`
)
{
    return I_0 * (1 - R - T) / (1 - T);  // same units as `I_0`
}
 
double nir_from_ppfd(
    double ppfd,                // micromol / m^2 / s
    double par_energy_content,  // J / micromol
    double par_energy_fraction  // dimensionless
)
{
    return ppfd * par_energy_content *
           (1 - par_energy_fraction) / par_energy_fraction;  // J / m^2 / s
}
 
double absorbed_shortwave(
    double incident_nir,            // J / m^2 / s
    double incident_ppfd,           // micromol / m^2 / s
    double par_energy_content,      // J / micromol
    double leaf_reflectance_par,    // dimensionless
    double leaf_transmittance_par,  // dimensionless
    double leaf_reflectance_nir,    // dimensionless
    double leaf_transmittance_nir   // dimensionless
)
{
    double const absorbed_par = thin_layer_absorption(
        leaf_reflectance_par,
        leaf_transmittance_par,
        incident_ppfd * par_energy_content);  // J / m^2 / s
 
    double const absorbed_nir = thin_layer_absorption(
        leaf_reflectance_nir,
        leaf_transmittance_nir,
        incident_nir);  // J / m^2 / s
 
    return absorbed_par + absorbed_nir;  // J / m^2 / s
}
 
double total_radiation(
    double Q_o,    // Light units such as `micromol / m^2 / s` or `J / m^2 / s`
    double k,      // dimensionless
    double alpha,  // dimensionless
    double ell     // dimensionless from m^2 leaf / m^2 ground
)
{
    return Q_o * exp(-k * sqrt(alpha) * ell);  // same units as `Q_ob`
}
 
double downscattered_radiation(
    double Q_ob,          // Light units such as `micromol / m^2 / s` or `J / m^2 / s`
    double k_direct,      // dimensionless
    double alpha_direct,  // dimensionless
    double ell            // dimensionless from m^2 leaf / m^2 ground
)
{
    return total_radiation(Q_ob, k_direct, alpha_direct, ell) -
           Q_ob * exp(-k_direct * ell);  // same units as `Q_ob`
}
 
double shaded_radiation(
    double Q_ob,       // Light units such as `micromol / m^2 / s` or `J / m^2 / s`
    double Q_od,       // same units as `Q_ob`
    double k_direct,   // dimensionless
    double k_diffuse,  // dimensionless
    double alpha,      // dimensionless
    double ell         // dimensionless from m^2 leaf / m^2 ground
)
{
    return total_radiation(Q_od, k_diffuse, alpha, ell) +
           downscattered_radiation(Q_ob, k_direct, alpha, ell);  // same units as `Q_ob`
}
 
Light_profile sunML(
    double ambient_ppfd_beam,       // micromol / (m^2 beam) / s
    double ambient_ppfd_diffuse,    // micromol / m^2 / s
    double chil,                    // dimensionless from m^2 / m^2
    double cosine_zenith_angle,     // dimensionless
    double heightf,                 // m^-1 from m^2 leaf / m^2 ground / m height
    double k_diffuse,               // dimensionless
    double lai,                     // dimensionless from m^2 / m^2
    double leaf_reflectance_nir,    // dimensionless
    double leaf_reflectance_par,    // dimensionless
    double leaf_transmittance_nir,  // dimensionless
    double leaf_transmittance_par,  // dimensionless
    double par_energy_content,      // J / micromol
    double par_energy_fraction,     // dimensionless
    int nlayers                     // dimensionless
)
{
    if (nlayers < 1 || nlayers > MAXLAY) {
        throw std::out_of_range("nlayers must be at least 1 but no more than MAXLAY.");
    }
 
    if (cosine_zenith_angle > 1 || cosine_zenith_angle < -1) {
        throw std::out_of_range("cosine_zenith_angle must be between -1 and 1.");
    }
 
    if (k_diffuse > 1 || k_diffuse < 0) {
        throw std::out_of_range("k_diffuse must be between 0 and 1.");
    }
 
    if (chil < 0) {
        throw std::out_of_range("chil must be non-negative.");
    }
 
    if (heightf <= 0) {
        throw std::out_of_range("heightf must greater than zero.");
    }
 
    // Calculate absorptivity from leaf reflectance and transmission
    double const absorptivity_nir = 1.0 - leaf_reflectance_nir - leaf_transmittance_nir; // dimensionless
    double const absorptivity_par = 1.0 - leaf_reflectance_par - leaf_transmittance_par; // dimensionless
 
    if (absorptivity_par > 1 || absorptivity_par < 0) {
        throw std::out_of_range("absorptivity_par must be between 0 and 1.");
    }
 
    if (absorptivity_nir > 1 || absorptivity_nir < 0) {
        throw std::out_of_range("absorptivity_nir must be between 0 and 1.");
    }
 
    // Calculate the leaf shape factor for an ellipsoidal leaf angle
    // distribution using the equation from page 251 of Campbell & Norman
    // (1998). We will use this value as `k_direct`, the canopy extinction
    // coefficient for direct photosynthetically active radiation throughout the
    // canopy. This quantity represents the ratio of horizontal area to total
    // area for leaves in the canopy and is therefore dimensionless from
    // (m^2 ground) / (m^2 leaf).
    double zenith_angle = acos(cosine_zenith_angle);  // radians
    double k0 = sqrt(pow(chil, 2) + pow(tan(zenith_angle), 2));
    double k1 = chil + 1.744 * pow((chil + 1.182), -0.733);
    double k_direct = k0 / k1;  // dimensionless
 
    double lai_per_layer = lai / nlayers;
 
    // Calculate the fraction of direct radiation that passes through the canopy
    // using Equation 15.1. Note that this is equivalent to the fraction of
    // ground area below the canopy that is exposed to direct sunlight. Note
    // that if the sun is at or below the horizon, no part of the soil is
    // sunlit; this corresponds to the case where cosine_zenith_angle is close
    // to or below zero.
    double canopy_direct_transmission_fraction =
        cosine_zenith_angle <= 1E-10 ? 0.0 : exp(-k_direct * lai);  // dimensionless
 
    // Calculate the ambient direct PPFD through a surface parallel to the ground
    const double ambient_ppfd_beam_ground = ambient_ppfd_beam * cosine_zenith_angle;  // micromol / (m^2 ground) / s
 
    // Calculate the ambient direct PPFD through a unit area of leaf surface
    double ambient_ppfd_beam_leaf = ambient_ppfd_beam_ground * k_direct;  // micromol / (m^2 leaf) / s
 
    // Calculate related NIR energy fluxes
    const double ambient_nir_beam = nir_from_ppfd(
        ambient_ppfd_beam, par_energy_content, par_energy_fraction);  // J / (m^2 beam) / s
 
    const double ambient_nir_beam_ground = nir_from_ppfd(
        ambient_ppfd_beam_ground, par_energy_content, par_energy_fraction);  // J / (m^2 ground) / s
 
    const double ambient_nir_diffuse = nir_from_ppfd(
        ambient_ppfd_diffuse, par_energy_content, par_energy_fraction);  // J / (m^2 ground) / s
 
    double ambient_nir_beam_leaf = nir_from_ppfd(
        ambient_ppfd_beam_leaf, par_energy_content, par_energy_fraction);  // J / (m^2 leaf) / s
 
    // Start to fill in the light profile values
    Light_profile light_profile;
    light_profile.canopy_direct_transmission_fraction = canopy_direct_transmission_fraction;
 
    // Fill in the layer-dependent light profile values
    for (int i = 0; i < nlayers; ++i) {
        // Get the cumulative LAI for this layer, which represents the total
        // leaf area above this layer
        const double cumulative_lai = lai_per_layer * (i + 0.5);
 
        // Calculate the PPFD incident on shaded leaves
        double shaded_ppfd = shaded_radiation(
            ambient_ppfd_beam_ground, ambient_ppfd_diffuse,
            k_direct, k_diffuse,
            absorptivity_par, cumulative_lai);  // micromol / m^2 / s
 
        // Calculate the NIR incident on shaded leaves
        double shaded_nir = shaded_radiation(
            ambient_nir_beam_ground, ambient_nir_diffuse,
            k_direct, k_diffuse,
            absorptivity_nir, cumulative_lai);  // J / m^2 / s
 
        // Calculate the fraction of sunlit and shaded leaves in this canopy
        // layer using Equation 15.22.
        double sunlit_fraction = exp(-k_direct * cumulative_lai);  // dimensionless
        double shaded_fraction = 1 - sunlit_fraction;              // dimensionless
 
        // For values of cosine_zenith_angle close to or less than 0, in place
        // of the calculations above, we want to use the limits of the above
        // expressions as cosine_zenith_angle approaches 0 from the right:
        if (cosine_zenith_angle <= 1E-10) {
            ambient_ppfd_beam_leaf = ambient_ppfd_beam / k1;
            ambient_nir_beam_leaf = ambient_nir_beam / k1;
            shaded_ppfd = ambient_ppfd_diffuse * exp(-k_diffuse * cumulative_lai);
            shaded_nir = ambient_nir_diffuse * exp(-k_diffuse * cumulative_lai);
            sunlit_fraction = 0;
            shaded_fraction = 1;
        }
 
        // Store values of incident PPFD
        light_profile.height[i] = (lai - cumulative_lai) / heightf;                    // m
        light_profile.shaded_fraction[i] = shaded_fraction;                            // dimensionless from m^2 / m^2
        light_profile.shaded_incident_ppfd[i] = shaded_ppfd;                           // micromol / (m^2 leaf) / s
        light_profile.shaded_incident_nir[i] = shaded_nir;                             // J / (m^2 leaf) / s
        light_profile.sunlit_fraction[i] = sunlit_fraction;                            // dimensionless from m^2 / m^2
        light_profile.sunlit_incident_ppfd[i] = ambient_ppfd_beam_leaf + shaded_ppfd;  // micromol / (m^2 leaf) / s
        light_profile.sunlit_incident_nir[i] = ambient_nir_beam_leaf + shaded_nir;     // J / (m^2 leaf) / s
 
        // Store values of absorbed PPFD
        light_profile.sunlit_absorbed_ppfd[i] =
            thin_layer_absorption(
                leaf_reflectance_par,
                leaf_transmittance_par,
                ambient_ppfd_beam_leaf + shaded_ppfd);  // micromol / m^2 / s
 
        light_profile.shaded_absorbed_ppfd[i] =
            thin_layer_absorption(
                leaf_reflectance_par,
                leaf_transmittance_par,
                shaded_ppfd);  // micromol / m^2 / s
 
        // Store values of absorbed solar energy (including PAR and NIR)
        light_profile.sunlit_absorbed_shortwave[i] =
            absorbed_shortwave(
                ambient_nir_beam_leaf + shaded_nir,
                ambient_ppfd_beam_leaf + shaded_ppfd,
                par_energy_content,
                leaf_reflectance_par,
                leaf_transmittance_par,
                leaf_reflectance_nir,
                leaf_transmittance_nir);  // J / (m^2 leaf) / s
 
        light_profile.shaded_absorbed_shortwave[i] =
            absorbed_shortwave(
                shaded_nir,
                shaded_ppfd,
                par_energy_content,
                leaf_reflectance_par,
                leaf_transmittance_par,
                leaf_reflectance_nir,
                leaf_transmittance_nir);  // J / (m^2 leaf) / s
    }
    return light_profile;
}