// root_system.C -- Root development and uptake.
// 
// Copyright 1996-2001 Per Abrahamsen and Søren Hansen
// Copyright 2000-2001 KVL.
//
// This file is part of Daisy.
// 
// Daisy is free software; you can redistribute it and/or modify
// it under the terms of the GNU Lesser Public License as published by
// the Free Software Foundation; either version 2.1 of the License, or
// (at your option) any later version.
// 
// Daisy is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU Lesser Public License for more details.
// 
// You should have received a copy of the GNU Lesser Public License
// along with Daisy; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

#define BUILD_DLL

#include "root_system.h"
#include "rootdens.h"
#include "librarian.h"
#include "geometry.h"
#include "soil_heat.h"
#include "soil_water.h"
#include "soil.h"
#include "chemical.h"
#include "chemistry.h"
#include "log.h"
#include "check.h"
#include "block_model.h"
#include "mathlib.h"
#include "librarian.h"
#include "treelog.h"
#include "frame.h"
#include "metalib.h"
#include <sstream>

double 
RootSystem::crown_potential () const
{ return h_x; }

double 
RootSystem::potential_water_uptake (const double h_x,
                                    const Geometry& geo,
                                    const Soil& soil,
                                    const SoilWater& soil_water,
                                    const double dt)
{
  const std::vector<double>& L = Density;
  std::vector<double>& S = H2OExtraction;
  const size_t size = std::min (S.size (), soil.size ());
  daisy_assert (L.size () >= size);
  const double area = M_PI * Rad * Rad; // [cm^2 R]

  for (size_t i = 0; i < size; i++)
    {
      if (L[i] <= 0.0 || soil_water.h (i) >= 0.0)
        {
          S[i] = 0.0;
          continue;
        }
      const double h = h_x - (1 + Rxylem) * geo.cell_z (i);
      daisy_assert (soil_water.Theta_ice (soil, i, h_wp) 
                    >= soil.Theta_res (i));
      const double max_uptake
        = std::max (0.0, (soil_water.Theta (i) 
                          - soil_water.Theta_ice (soil, i, h_wp)) / dt);
      const double uptake
        = bound (0.0, 
                 (2 * M_PI * L[i]
                  * (soil_water.Theta_ice (soil, i, h) 
                     / soil_water.Theta_ice (soil, i, 0.0))
                  * (soil.M (i, soil_water.h (i)) - soil.M (i, h))
                  / (- 0.5 * log (area * L[i]))),
                 max_uptake);
      daisy_assert (soil_water.h (i) > h_wp || iszero (uptake));
      daisy_assert (L[i] >= 0.0);
      daisy_assert (soil_water.Theta_ice (soil, i, h) > 0.0);
      daisy_assert (soil_water.Theta_ice (soil, i, 0.0) > 0.0);
      daisy_assert (soil.M (i, soil_water.h (i)) >= 0.0);
      daisy_assert (soil.M (i, h) >= 0.0);
      daisy_assert (area * L[i] > 0.0);
      daisy_assert (std::isnormal (- 0.5 * log (area * L[i])));
      daisy_assert (uptake >= 0.0);
      S[i] = uptake;
    }
  return geo.total_surface (S) * 10.0 /* [mm/cm] */;
}

double
RootSystem::water_uptake (const Units& units, double Ept_,
                          const Geometry& geo,
                          const Soil& soil,
                          const SoilWater& soil_water,
                          const double EvapInterception,
                          const double dt,
                          Treelog& msg)
{
  daisy_assert (EvapInterception >= 0);
  if (Ept_ < 0)
    {
      Treelog::Open nest (msg, "RootSystem water uptake");
      std::ostringstream tmp;
      tmp << "BUG: Negative EPT (" << Ept_ << ")";
      msg.error (tmp.str ());
      Ept_ = 0.0;
    }
  Ept = Ept_;

  static const double min_step = 1.0;
  double total = potential_water_uptake (h_x, geo, soil, soil_water, dt);
  double step = min_step;

  while (total < Ept && h_x > h_wp)
    {
      const double h_next = std::max (h_x - step, h_wp);
      const double next = potential_water_uptake (h_next, geo, 
                                                  soil, soil_water, dt);

      if (next < total)
        {
          // We are past the top of the curve.
          if (step <= min_step)
            // We cannot go any closer to the top, skip it.
            {
              h_x = h_wp;
              total = potential_water_uptake (h_x, geo, soil, soil_water, dt);
              break;
            }
          else
            // Try again a little close.
            {
              step /= 2;
              continue;
            }
        }
      total = next;
      h_x = h_next;
      step *= 2;
    }
  if (h_x < h_wp)
    h_x = h_wp;

  step = min_step;
  daisy_assert (h_x < 0.001);
  while (total > Ept && h_x < 0.0)
    {
      daisy_assert (h_x < 0.001);
      const double h_next = std::min (h_x + step, 0.0);
      const double next = potential_water_uptake (h_next, geo, 
                                                  soil, soil_water, dt);

      if (next < Ept)
        {
          // We went too far.
          if (step <= min_step)
            {
              // We can't get any closer.
              daisy_assert (next <= total);
              if (next >= Ept)
                {
                  // total = next;
                  h_x = h_next;
                }
              else

                break;
            }
          else
            // Try again a little closer.
            {
              step /= 2;
              continue;
            }
        }
      total = next;
      h_x = h_next;
      step *= 2;
    }

  // We need this to make sure H2OExtraction corresponds to 'h_x'.
  const double total2 
    = potential_water_uptake (h_x, geo, soil, soil_water, dt);
  daisy_assert (approximate (total, total2));
  daisy_assert (h_x >= h_wp);

  if (total > Ept)
    {
      daisy_assert (h_x < 0.001);
      daisy_assert (total > 0);
      const double factor = Ept / total;
      for (size_t i = 0; i < soil.size (); i++)
        H2OExtraction[i] *= factor;
      total = Ept;
    }
  H2OUpt = total;

  // Update water stress factor
  if (Ept < 0.010)
    water_stress = 0.0;
  else
    water_stress = 1.0 - (total + EvapInterception) / (Ept + EvapInterception);

  // ABA production.
  ABAprod->production (units, geo, soil_water, H2OExtraction, Density,
                       ABAExtraction, msg);
  const double mm_per_cm = 10.0; // [mm/cm]
  if (H2OUpt > 0.0)
    // [g/cm^3 W] = [g/cm^2 A] * [mm/cm] / [mm/h]
    ABAConc = geo.total_surface (ABAExtraction) * mm_per_cm / H2OUpt;
  else
    { /* Use old value */ }
  daisy_assert (std::isfinite (ABAConc));

  // Result.
  return H2OUpt;
}

double
RootSystem::solute_uptake (const Geometry& geo, const Soil& soil,
                           const SoilWater& soil_water,
                           Chemical& solute,
                           double PotNUpt,
                           std::vector<double>& uptake,
                           const double I_max,      // [g/cm R/h]
                           const double C_root_min, // [g/cm^3 W]
                           const double dt)
{
  if (PotNUpt <= 0.0)
    {
      fill (uptake.begin (), uptake.end (), 0.0);
      return 0.0;
    }

  daisy_assert (PotNUpt > 0.0);
  PotNUpt /= 1.0e4;             // g/m^2/h -> g/cm^2/h
  PotNUpt *= geo.surface_area (); // [g/h]
  const int size = soil.size ();
  // I: Uptake per root length.
  // I = I_zero - B_zero * C_root
  std::vector<double> I_zero (size, 0.0); // [g/cm R/h]
  std::vector<double> B_zero (size, 0.0); // [cm^3 W/cm R/h] 

  double U_zero = 0.0;          // [g/h] Total uptake a C_root_min
  double B = 0.0;               // [cm^3 W/h]

  for (int i = 0; i < size; i++)
    {
      const double C_l = solute.C_secondary (i);     // [g/cm^3 W]
      const double Theta = soil_water.Theta_old (i); // [cm^3 W/cm^3 S]
      const double L = Density[i];                   // [cm R/cm^3 S]
      if (L > 0 && soil_water.h (i) <= 0.0)
        {
          // [cm^3 W/cm R/h] = [cm^3 W/cm^3 S/h] / [cm R/cm^3 S]
          const double q_r = H2OExtraction[i] / L;
          // [cm^2/h] = [cm^2/h] * [cm^3/cm^3] * [cm^3 W/cm^3 S]
          const double D = solute.diffusion_coefficient ()
            * soil.tortuosity_factor (i, Theta)
            * Theta;
          // [] = [cm^3 W/cm R/h] / [cm^2/h]
          const double alpha = q_r / ( 2 * M_PI * D);
          // [] = [] / ([cm] * sqrt ([cm^-2]))
          const double beta = 1.0 / (Rad * sqrt (M_PI * L));
          // [] = [] * []
          const double beta_squared = beta * beta;
          if (alpha < 1e-10)
            {
              // [cm^3 W/cm R/h] = [cm^2/h] / ([] * log ([]/[]) - [])
              B_zero[i] = 4.0 * M_PI * D
                / (beta_squared * log (beta_squared) / (beta_squared - 1.0)
                   - 1.0);
              // [g/cm R/h] = [cm^3 W/cm R/h] * [g/cm^3 W]
              I_zero[i] = B_zero[i] * C_l;
            }
          else
            { 
              const double divisor // []
                = ((beta_squared - 1.0) * (1.0 - 0.5 * alpha)
                   - (pow (beta, 2.0 - alpha) - 1.0));
              
              if (std::isnormal (divisor))
                {
                  // [cm^3 W/cm R/h] = [cm^3 W/cm R/h]
                  B_zero[i] = q_r * (pow (beta, 2.0 - alpha) - 1.0)
                    / divisor;
                  // [g/cm R/h] = [cm^3 W/cm R/h] * [g/cm^3 W]
                  I_zero[i] = q_r * (beta_squared - 1.0) * (1.0 - 0.5 * alpha)
                    * C_l / divisor;
                }
              else              
                {
                  daisy_assert (approximate (alpha, 2.0));
                  const double div2 // []
                    = ((beta_squared - 1.0) - log (beta_squared));
                  daisy_assert (std::isnormal (div2));
                  // [cm^3 W/cm R/h] = [cm^3 W/cm R/h]
                  B_zero[i] = q_r * log (beta_squared) / div2;
                  // [g/cm R/h] = [cm^3 W/cm R/h] * [g/cm^3 W]
                  I_zero[i] = q_r * (beta_squared - 1.0) * C_l / div2;
                }
            }
          daisy_assert (std::isfinite (I_zero[i]));
          daisy_assert (std::isfinite (B_zero[i]));
          // [cm^3 W/h] = [cm R/cm^3 S] * [cm^3 S] * [cm^3 W/cm R/h]
          B += L * geo.cell_volume (i) * B_zero[i];
          // [g/h] = [cm R/cm^3 S] * [cm^3 S]
          //       * ([g/cm R/h] - [cm^3 W/cm R/h] * [g/cm^3 W], [g/cm R/h])
          U_zero += L * geo.cell_volume (i) 
            * bound (0.0, I_zero[i] - B_zero[i] * C_root_min, I_max);
        }
    }
  double C_root = C_root_min;   // [g/cm^3 W]
  if (U_zero > PotNUpt)
    // [g/cm^3 W] = max (([g/h] - [g/h]) / [cm^3 W/h], [g/cm^3 W])
    C_root = std::max ((U_zero - PotNUpt) / B, C_root_min);

  for (int i = 0; i < size; i++)
    {
      const double L = Density[i]; // [cm R/cm^3 S]
      if (L > 0 && soil_water.h (i) <= 0.0)
        // [g/cm^3 S/h]
        //  = [cm R/cm^3 S] * ([g/cm R/h] - [cm^3 W/cm R/h] * [g/cm^3 W])
        uptake[i] = std::max (0.0, L * (std::min (I_zero[i], I_max)
                                        - B_zero[i] * C_root));
      else
        uptake[i] = 0.0;
      daisy_assert (uptake[i] >= 0.0);
    }
  solute.add_to_root_sink (uptake);

  // gN/cm³/h -> gN/m²/h
  return geo.total_surface (uptake) * 1.0e4;
}

double
RootSystem::nitrogen_uptake (const Geometry& geo, const Soil& soil,
                             const SoilWater& soil_water,
                             Chemistry& chemistry,
                             const double NH4_root_min,
                             const double NO3_root_min,
                             const double PotNUpt, double dt)
{
  if (PotNUpt <= 0.0)           
    // No uptake needed.
    {
      NH4Upt = NO3Upt = 0.0;
      fill (NH4Extraction.begin (), NH4Extraction.end (), 0.0);
      fill (NO3Extraction.begin (), NO3Extraction.end (), 0.0);
    }
  else if (chemistry.know (Chemical::NH4 ()) 
           && chemistry.know (Chemical::NO3 ()))
    // Normlal uptake.
    {
      Chemical& soil_NH4 = chemistry.find (Chemical::NH4 ());
      Chemical& soil_NO3 = chemistry.find (Chemical::NO3 ());

      NH4Upt = solute_uptake (geo, soil, soil_water, soil_NH4, 
                              PotNUpt, NH4Extraction,
                              MxNH4Up, NH4_root_min, dt);
      NO3Upt = solute_uptake (geo, soil, soil_water, soil_NO3, 
                              PotNUpt - NH4Upt, NO3Extraction, 
                              MxNO3Up, NO3_root_min, dt);
    }
  else 
    // If we don't track inorganic N, assume we have enough.
    {
      NH4Upt = PotNUpt;
      NO3Upt = 0.0;
    }
  daisy_assert (NH4Upt >= 0.0);
  daisy_assert (NO3Upt >= 0.0);

  return NH4Upt + NO3Upt;
}

void
RootSystem::tick_dynamic (const double T, const double day_fraction,
                          SoilWater& soil_water, const double dt)
{
  // Update soil water sink term.
  soil_water.root_uptake (H2OExtraction);

  // Accumulated water stress.
  water_stress_days += water_stress * day_fraction;

  // Keep track of daily soil temperature.
  partial_soil_temperature += T * dt;
  partial_day += dt;
  if (partial_day >= 24.0)
    {
      soil_temperature = partial_soil_temperature / partial_day;
      partial_soil_temperature = 0.0;
      partial_day = 0.0;
    }

  // Clear nitrogen.
  NH4Upt = NO3Upt =0.0;
}

void
RootSystem::tick_daily (const Geometry& geo, const Soil& soil, 
                        const double WRoot, const bool root_growth,
                        const double DS, Treelog& msg)
{
  const double SoilLimit = -soil.MaxRootingHeight ();
  // Penetration.
  if (root_growth)
    {
      const double clay = geo.content_height (soil, &Soil::clay, -Depth);
      double clay_fac = PenClayFac (clay);
      double dp = PenPar1 * clay_fac * std::max (0.0, soil_temperature - PenPar2);
      PotRtDpt = std::min (PotRtDpt + dp, MaxPen);
      /*max depth determined by crop*/
      Depth = std::min (Depth + dp, MaxPen);
      PotRtDpt = std::max (PotRtDpt, Depth);
      /*max depth determined by crop*/
      Depth = std::min (Depth, SoilLimit); /*or by soil conditions*/
    }
  set_density (geo, SoilLimit, WRoot, DS, msg);
}

void
RootSystem::set_density (const Geometry& geo, const double SoilLimit,
                         const double WRoot, const double DS, Treelog& msg)
{ rootdens->set_density (geo, SoilLimit, PotRtDpt, 
                         PotRtDpt * (MaxWidth / MaxPen),
                         WRoot, DS, Density, msg); }

void
RootSystem::full_grown (const Geometry& geo, 
                        const double max_rooting_depth,
                        const double WRoot, Treelog& msg)
{
  PotRtDpt = MaxPen;
  Depth = std::min (MaxPen, -max_rooting_depth);
  set_density (geo, -max_rooting_depth, WRoot, 1.0, msg);
}

void
RootSystem::output (Log& log) const
{
  output_object (rootdens, "rootdens", log);
  output_derived (ABAprod, "ABAprod", log);
  output_variable (PotRtDpt, log);
  output_variable (Depth, log);
  output_variable (Density, log);
  output_variable (H2OExtraction, log);
  output_variable (NH4Extraction, log);
  output_variable (NO3Extraction, log);
  output_variable (ABAExtraction, log);
  output_variable (ABAConc, log);
  output_variable (h_x, log);
  output_variable (partial_soil_temperature, log);
  output_variable (partial_day, log);
  output_variable (soil_temperature, log);
  output_variable (water_stress, log);
  output_variable (water_stress_days, log);
  output_variable (production_stress, log);
  output_variable (Ept, log);
  output_variable (H2OUpt, log);
  output_variable (NH4Upt, log);
  output_variable (NO3Upt, log);
}

void
RootSystem::initialize (const Metalib& metalib, 
                        const Geometry& geo, const double row_width, 
                        const double row_pos, Treelog& msg)
{
  const bool is_row_crop = row_width > 0.0;
  if (rootdens.get ())
    /* We already has a root density model. */;
  else if (is_row_crop)
    rootdens = Rootdens::create_row (metalib, msg, row_width, row_pos);
  else
    rootdens = Rootdens::create_uniform (metalib, msg);

  rootdens->initialize (geo, row_width, row_pos, msg);
  initialize (metalib.units (), geo, msg);
}

void
RootSystem::initialize (const Units& units, const Geometry& geo, Treelog& msg)
{
  const size_t cell_size = geo.cell_size ();
  ABAprod->initialize (units, msg);
  while (Density.size () < cell_size)
    Density.push_back (0.0);
  while (H2OExtraction.size () < cell_size)
    H2OExtraction.push_back (0.0);
  while (NH4Extraction.size () < cell_size)
    NH4Extraction.push_back (0.0);
  while (NO3Extraction.size () < cell_size)
    NO3Extraction.push_back (0.0);
  while (ABAExtraction.size () < cell_size)
    ABAExtraction.push_back (0.0);
}

bool
RootSystem::check (const Units& units, Treelog& msg) const
{
  bool ok = true;
  if (!ABAprod->check (units, msg))
    ok = false;
  return ok;
}

void 
RootSystem::load_syntax (Frame& frame)
{

  frame.declare_object ("rootdens", Rootdens::component, 
                        Attribute::OptionalConst, Attribute::Singleton,
                        "Root density model.");

  frame.declare_object ("ABAprod", ABAProd::component, "ABA production model.");
  frame.set ("ABAprod", "none");

  frame.declare ("DptEmr", "cm", Check::non_negative (), Attribute::Const,
                 "Penetration at emergence.");
  frame.set ("DptEmr", 10.0);
  frame.declare ("PenPar1", "cm/dg C/d", Check::non_negative (), Attribute::Const,
                 "Penetration rate parameter, coefficient.");
  frame.set ("PenPar1", 0.25);
  frame.declare ("PenPar2", "dg C", Check::none (), Attribute::Const,
                 "Penetration rate parameter, threshold.");
  frame.set ("PenPar2", 4.0);
  frame.declare ("PenClayFac", Attribute::Fraction (), Attribute::None (),
                 Check::non_negative (), Attribute::Const, 
                 "Clay dependent factor to multiply 'PenPar1' with.");
  PLF clay;
  clay.add (0.0, 1.0);
  clay.add (1.0, 1.0);
  frame.set ("PenClayFac", clay);
  frame.declare ("MaxPen", "cm", Check::positive (), Attribute::Const,
                 "Maximum penetration depth.");
  frame.set ("MaxPen", 100.0);
  frame.declare ("MaxWidth", "cm", Check::positive (), Attribute::OptionalConst,
                 "Maximum horizontal distance of roots from plant.");
  frame.declare ("Rad", "cm", Check::positive (), Attribute::Const,
                 "Root radius.");
  frame.set ("Rad", 0.005);
  frame.declare ("h_wp", "cm", Check::none (), Attribute::Const,
                 "Matrix potential at wilting point.");
  frame.set ("h_wp",-15000.0);
  frame.declare ("MxNH4Up", "g/cm/h", Check::non_negative (), Attribute::Const,
                 "Maximum NH4 uptake per unit root length.");
  frame.set ("MxNH4Up", 2.5e-7);
  frame.declare ("MxNO3Up", "g/cm/h", Check::non_negative (), Attribute::Const,
                 "Maximum NO3 uptake per unit root length.");
  frame.set ("MxNO3Up", 2.5e-7);
  frame.declare ("Rxylem", Attribute::None (), Check::non_negative (), Attribute::Const,
                 "Transport resistence in xyleme.");
  frame.set ("Rxylem", 10.0);

  frame.declare ("PotRtDpt", "cm", Check::non_negative (), Attribute::OptionalState,
                 "Potential root penetration depth.");
  frame.declare ("Depth", "cm", Check::non_negative (), Attribute::OptionalState,
                 "Rooting Depth.");
  frame.declare ("Density", "cm/cm^3", Check::non_negative (),
                 Attribute::LogOnly, Attribute::SoilCells,
                 "Root density in soil layers.");
  frame.declare ("H2OExtraction", "cm^3/cm^3/h", Check::non_negative (), 
                 Attribute::LogOnly, Attribute::SoilCells,
                 "Extraction of H2O in soil layers.");
  frame.declare ("NH4Extraction", "g N/cm^3/h", Check::non_negative (), 
                 Attribute::LogOnly, Attribute::SoilCells,
                 "Extraction of NH4-N in soil layers.");
  frame.declare ("NO3Extraction", "g N/cm^3/h", Check::non_negative (), 
                 Attribute::LogOnly, Attribute::SoilCells,
                 "Extraction of NO3-N in soil layers.");
  frame.declare ("ABAExtraction", "g/cm^3/h", Check::non_negative (), 
                 Attribute::LogOnly, Attribute::SoilCells,
                 "Extraction of ABA in soil layers.");
  frame.declare ("ABAConc", "g/cm^3", Check::non_negative (), 
                 Attribute::State, "ABA concentration in water uptake.");
  frame.set ("ABAConc", 0.0);
  frame.declare ("h_x", "cm", Check::none (), Attribute::State,
                 "Root extraction at surface.");
  frame.set ("h_x", 0.0);
  frame.declare ("partial_soil_temperature", "dg C h", Attribute::State,
                 "Soil temperature hours this day, so far.");
  frame.set ("partial_soil_temperature", 0.0);
  frame.declare ("partial_day", "h", Attribute::State,
                 "Hours we have accumulated soil temperature this day.");
  frame.set ("partial_day", 0.0);
  frame.declare ("soil_temperature", "dg C", Attribute::State,
                 "Average soil temperature yesterday.");
  frame.set ("soil_temperature", 0.0);
  frame.declare ("water_stress", Attribute::None (), Check::fraction (),
                 Attribute::LogOnly,
                 "Fraction of requested water we didn't get.");
  frame.declare ("water_stress_days", "d", Check::non_negative (),
                 Attribute::State,
                 "Number of days production has halted due to water stress.\n\
This is the sum of water stress for each hour, multiplied with the\n\
fraction of the radition of that day that was received that hour.");
  frame.set ("water_stress_days", 0.0);
  frame.declare ("production_stress", Attribute::None (), Check::fraction (), 
                 Attribute::LogOnly,
                 "SVAT induced stress, or -1 if not applicable.");
  frame.declare ("Ept", "mm/h", Check::none (), Attribute::LogOnly,
                 "Potential transpiration.");
  frame.declare ("H2OUpt", "mm/h", Check::non_negative (), Attribute::LogOnly,
                 "H2O uptake.");
  frame.declare ("NH4Upt", "g N/m^2/h", Check::non_negative (), Attribute::LogOnly,
                 "NH4-N uptake.");
  frame.declare ("NO3Upt", "g N/m^2/h", Check::non_negative (), Attribute::LogOnly,
                 "NO3-N uptake.");
}

static double
get_PotRtDpt (const Block& al)
{
  if (al.check ("PotRtDpt"))
    return al.number ("PotRtDpt");
  if (al.check ("Depth"))
    return al.number ("Depth");
  return al.number ("DptEmr");
}

RootSystem::RootSystem (const Block& al)
  : rootdens (al.check ("rootdens") 
              ? Librarian::build_item<Rootdens> (al, "rootdens")
              : NULL),
    ABAprod (Librarian::build_item<ABAProd> (al, "ABAprod")),
    PenPar1 (al.number ("PenPar1")),
    PenPar2 (al.number ("PenPar2")),
    PenClayFac (al.plf ("PenClayFac")),
    MaxPen (al.number ("MaxPen")),
    MaxWidth (al.number ("MaxWidth", MaxPen)),
    Rad (al.number ("Rad")),
    h_wp (al.number ("h_wp")),
    MxNH4Up (al.number ("MxNH4Up")),
    MxNO3Up (al.number ("MxNO3Up")),
    Rxylem (al.number ("Rxylem")),
    PotRtDpt (get_PotRtDpt (al)),
    Depth (al.check ("Depth") ? al.number ("Depth") : al.number ("DptEmr")),
    ABAConc (al.number ("ABAConc")),
    h_x (al.number ("h_x")),
    partial_soil_temperature (al.number ("partial_soil_temperature")),
    partial_day (al.number ("partial_day")),
    soil_temperature (al.number ("soil_temperature")),
    water_stress (0.0),
    water_stress_days (al.number ("water_stress_days")),
    production_stress (-1.0),
    Ept (0.0),
    H2OUpt (0.0),
    NH4Upt (0.0),
    NO3Upt (0.0)
{ }

RootSystem::~RootSystem ()
{ }

static DeclareSubmodel 
root_system_submodel (RootSystem::load_syntax, "RootSystem", "\
Standard root system model.");
