bool unique_validate (const std::vector<T>& list, Treelog& msg)
{
  bool ok = true;
  std::map<T, size_t> found;
  for (size_t i = 0; i < list.size (); i++)
    {
      T v = list[i];
      typename std::map<T, size_t>::const_iterator f = found.find (v);
      if (f != found.end ())
	{
	  std::ostringstream tmp;
	  tmp << "Entry " << (*f).second << " and " << i 
	      << " are both '" << v << "'";
	  msg.error (tmp.str ());
          ok = false;
	}
      found[v] = i;
    }
  return ok;
}

void
SoilWater::drain (const std::vector<double>& v, Treelog& msg)
{
  forward_sink (v, v);

  daisy_assert (S_sum_.size () == v.size ());
  daisy_assert (S_drain_.size () == v.size ());
  daisy_assert (S_soil_drain_.size () == v.size ());
  for (unsigned i = 0; i < v.size (); i++)
    {
      if (v[i] < -1e-8)
        {
          std::ostringstream tmp;
          tmp << "draining " << v[i] << " [h^-1] from cell " << i;
          msg.debug (tmp.str ());
        }
      S_sum_[i] += v[i];
      S_drain_[i] += v[i];
      S_soil_drain_[i] += v[i];
    }
}

void 
ChemistryMulti::deposit (const symbol chem, const double flux, Treelog& msg)
{
  bool found = false;

  for (size_t c = 0; c < combine.size (); c++)
    if (combine[c]->know (chem))
      {
	if (found)
	  msg.error ("Duplicate chemical '" + chem + "' detected");

	Chemical& chemical = combine[c]->find (chem);
        chemical.deposit (flux);
	found = true;
      }
  
  if (found)
    return;

  check_ignore (chem, msg);
}

bool 
Frame::check (const Metalib& metalib, const Frame& frame,
              const symbol key, Treelog& msg) const
{ 
  if (lookup (key) == Attribute::Error)
    {
      msg.error ("'" + key + "': not defined");
      return false;
    }

  const Type& type = impl->find_type (frame, key);
  bool ok = true;

  if (!impl->check (metalib, frame, type, key, msg))
    ok = false;
  else if (parent ()
      && parent ()->lookup (key) != Attribute::Error
      && !parent ()->check (metalib, frame, key, msg))
    ok = false;

  return ok;
}

void
SoilWater::initialize (const FrameSubmodel& al, const Geometry& geo,
                       const Soil& soil, const SoilHeat& soil_heat,
                       const Groundwater& groundwater, Treelog& msg)
{
  Treelog::Open nest (msg, "SoilWater");

  const size_t cell_size = geo.cell_size ();
  const size_t edge_size = geo.edge_size ();

  // Ice must be first.
  if (al.check ("X_ice"))
    {
      X_ice_ = al.number_sequence ("X_ice");
      if (X_ice_.size () == 0)
        X_ice_.push_back (0.0);
      while (X_ice_.size () < cell_size)
        X_ice_.push_back (X_ice_[X_ice_.size () - 1]);
    }
  else 
    X_ice_.insert (X_ice_.begin (), cell_size, 0.0);

  if (al.check ("X_ice_buffer"))
    {
      X_ice_buffer_ = al.number_sequence ("X_ice_buffer");
      if (X_ice_buffer_.size () == 0)
        X_ice_buffer_.push_back (0.0);
      while (X_ice_buffer_.size () < cell_size)
        X_ice_buffer_.push_back (X_ice_buffer_[X_ice_buffer_.size () - 1]);
    }
  else 
    X_ice_buffer_.insert (X_ice_buffer_.begin (), cell_size, 0.0);

  for (size_t i = 0; i < cell_size; i++)
    {
      const double Theta_sat = soil.Theta (i, 0.0, 0.0);
      daisy_assert (Theta_sat >= X_ice_[i]);
      h_ice_.push_back (soil.h (i, Theta_sat - X_ice_[i]));
    }
  daisy_assert (h_ice_.size () == cell_size);

  geo.initialize_layer (Theta_, al, "Theta", msg);
  geo.initialize_layer (h_, al, "h", msg);

  for (size_t i = 0; i < Theta_.size () && i < h_.size (); i++)
    {
      const double Theta_h = soil.Theta (i, h_[i], h_ice (i));
      if (!approximate (Theta_[i], Theta_h))
        {
          std::ostringstream tmp;
          tmp << "Theta[" << i << "] (" << Theta_[i] << ") != Theta (" 
              << h_[i] << ") (" << Theta_h << ")";
          msg.error (tmp.str ());
        }
      Theta_[i] = Theta_h;
    }
  if (Theta_.size () > 0)
    {
      while (Theta_.size () < cell_size)
        Theta_.push_back (Theta_[Theta_.size () - 1]);
      if (h_.size () == 0)
        for (size_t i = 0; i < cell_size; i++)
          h_.push_back (soil.h (i, Theta_[i]));
    }
  if (h_.size () > 0)
    {
      while (h_.size () < cell_size)
        h_.push_back (h_[h_.size () - 1]);
      if (Theta_.size () == 0)
        for (size_t i = 0; i < cell_size; i++)
          Theta_.push_back (soil.Theta (i, h_[i], h_ice (i)));
    }
  daisy_assert (h_.size () == Theta_.size ());

  // Groundwater based pressure.
  if (h_.size () == 0)
    {
      if (groundwater.table () > 0.0)
        {
          const double h_pF2 = -100.0; // pF 2.0;
          for (size_t i = 0; i < cell_size; i++)
            {
              h_.push_back (h_pF2);
              Theta_.push_back (soil.Theta (i, h_pF2, h_ice_[i]));
            }
        }
      else
        {
          const double table = groundwater.table ();
          
          for (size_t i = 0; i < cell_size; i++)
            {
              h_.push_back (std::max (-100.0, table - geo.cell_z (i)));
              Theta_.push_back (soil.Theta (i, h_[i], h_ice_[i]));
            }
        }
    }
  daisy_assert (h_.size () == cell_size);

  // Sources.
//.........这里部分代码省略.........

void
SoilWater::tick_after (const Geometry& geo,
                       const Soil& soil, const SoilHeat& soil_heat, 
                       const bool initial, 
                       Treelog& msg)
{
  TREELOG_SUBMODEL (msg, "SoilWater");

  // We need old K for primary/secondary flux division.
  std::vector<double> K_old = K_cell_;

  // Update cells.
  const size_t cell_size = geo.cell_size ();
  daisy_assert (K_cell_.size () == cell_size);
  daisy_assert (Cw2_.size () == cell_size);
  daisy_assert (h_.size () == cell_size);
  daisy_assert (h_ice_.size () == cell_size);
  daisy_assert (K_old.size () == cell_size);
  daisy_assert (Theta_.size () == cell_size);
  daisy_assert (Theta_primary_.size () == cell_size);
  daisy_assert (Theta_secondary_.size () == cell_size);
  daisy_assert (Theta_tertiary_.size () == cell_size);

  double z_low = geo.top ();
  table_low = NAN;
  double z_high = geo.bottom ();
  table_high = NAN;

  for (size_t c = 0; c < cell_size; c++)
    {
      // Groundwater table.
      const double z = geo.cell_top (c);
      const double h = h_[c];
      const double table = z + h;
      if (h < 0)
        {
          if (approximate (z, z_low))
            {
              if (!std::isnormal (table_low)
                  || table < table_low)
                table_low = table;
            }
          else if (z < z_low)
            table_low = table;
        }
      else if (approximate (z, z_high))
        {
          if (!std::isnormal (table_high)
              || table > table_high)
            table_high = table;
        }
      else if (z > z_high)
        table_high = table;

      // Conductivity.
      K_cell_[c] = soil.K (c, h_[c], h_ice_[c], soil_heat.T (c));
      
      // Specific water capacity.
      Cw2_[c] = soil.Cw2 (c, h_[c]);

      // Primary and secondary water.
      if (Theta_[c] <= 0.0)
        {
          std::ostringstream tmp;
          tmp << "Theta[" << c << "] = " << Theta_[c];
          daisy_bug (tmp.str ());
          Theta_[c] = 1e-9;
        }
      const double h_lim = soil.h_secondary (c);
      if (h_lim >= 0.0 || h_[c] <= h_lim)
        {
          // No active secondary domain.
          Theta_primary_[c] = Theta_[c];
          Theta_secondary_[c] = 0.0;
        }
      else 
        {
          // Secondary domain activated.
          const double Theta_lim = soil.Theta (c, h_lim, h_ice_[c]);
          daisy_assert (Theta_lim > 0.0);
          if (Theta_[c] >= Theta_lim)
            {
              Theta_primary_[c] = Theta_lim;
              Theta_secondary_[c] = Theta_[c] - Theta_lim;
            }
          else
            { 
              std::ostringstream tmp;
              tmp << "h[" << c << "] = " << h_[c] 
                  << "; Theta[" << c << "] = " << Theta_[c] 
                  << "\nh_lim = " << h_lim << "; Theta_lim = " << Theta_lim
                  << "\nStrenge h > h_lim, yet Theta <= Theta_lim";
              msg.bug (tmp.str ());
              Theta_primary_[c] = Theta_[c];
              Theta_secondary_[c] = 0.0;
            }
        }
    }

  if (!std::isnormal (table_high))
//.........这里部分代码省略.........

void 
Movement1D::tick_water (const Geometry1D& geo,
                        const Soil& soil, const SoilHeat& soil_heat, 
                        Surface& surface, Groundwater& groundwater,
                        const std::vector<double>& S,
                        std::vector<double>& h_old,
                        const std::vector<double>& Theta_old,
                        const std::vector<double>& h_ice,
                        std::vector<double>& h,
                        std::vector<double>& Theta,
                        std::vector<double>& q,
                        std::vector<double>& q_p,
                        const double dt, 
                        Treelog& msg)
{
  const size_t top_edge = 0U;
  const size_t bottom_edge = geo.edge_size () - 1U;

  // Limit for groundwater table.
  size_t last  = soil.size () - 1;

  // Limit for ridging.
  const size_t first = (surface.top_type (geo, 0U) == Surface::soil)
    ? surface.last_cell (geo, 0U) 
    : 0U;
  // Calculate matrix flow next.

  for (size_t m = 0; m < matrix_water.size (); m++)
    {
      water_attempt (m);
      Treelog::Open nest (msg, matrix_water[m]->name);
      try
        {
          matrix_water[m]->tick (msg, geo, soil, soil_heat,
                                 first, surface, top_edge, 
                                 last, groundwater, bottom_edge,
                                 S, h_old, Theta_old, h_ice, h, Theta, 0U, q, 
                                 dt);

          for (size_t i = last + 2; i <= soil.size (); i++)
            {
              q[i] = q[i-1];
              q_p[i] = q_p[i-1];
            }

          // Update surface and groundwater reservoirs.
          surface.accept_top (q[0] * dt, geo, 0U, dt, msg);
          surface.update_pond_average (geo);
          const double q_down = q[soil.size ()] + q_p[soil.size ()];
          groundwater.accept_bottom (q_down * dt, geo, soil.size ());
          if (m > 0)
            msg.debug ("Reserve model succeeded");
          return;
        }
      catch (const char* error)
        {
          msg.debug (std::string ("UZ problem: ") + error);
        }
      catch (const std::string& error)
        {
          msg.debug (std::string ("UZ trouble: ") + error);
        }
      
      water_failure (m);
    }
  throw "Water matrix transport failed"; 
}

 void doIt (Daisy&, const Scope&, Treelog& out)
 { 
   out.warning (message.name ());
 }

 void doIt (Daisy& daisy, const Scope&, Treelog& out)
   { 
     out.message ("Ridging");      
     daisy.field ().ridge (ridge); 
   }

void
MovementSolute::solute (const Soil& soil, const SoilWater& soil_water,
                        const double J_above, Chemical& chemical, 
                        const double dt,
                        const Scope& scope, Treelog& msg)
{
  daisy_assert (std::isfinite (J_above));
  const size_t cell_size = geometry ().cell_size ();
  const size_t edge_size = geometry ().edge_size ();

  // Source term transfered from secondary to primary domain.
  std::vector<double> S_extra (cell_size, 0.0);

  // Divide top solute flux according to water.
  std::map<size_t, double> J_tertiary;
  std::map<size_t, double> J_secondary; 
  std::map<size_t, double> J_primary;

  if (J_above > 0.0)
    // Outgoing, divide according to content in primary domain only.
    divide_top_outgoing (geometry (), chemical, J_above, 
                         J_primary, J_secondary, J_tertiary);
  else if (J_above < 0.0)
    // Incomming, divide according to all incomming water.
    divide_top_incomming (geometry (), soil_water, J_above, 
                          J_primary, J_secondary, J_tertiary);
  else
    // No flux.
    zero_top (geometry (), J_primary, J_secondary, J_tertiary);

  // Check result.
  {
    const std::vector<size_t>& edge_above 
      = geometry ().cell_edges (Geometry::cell_above);
    const size_t edge_above_size = edge_above.size ();
    double J_sum = 0.0;
    for (size_t i = 0; i < edge_above_size; i++)
      {
        const size_t edge = edge_above[i];
        const double in_sign 
          = geometry ().cell_is_internal (geometry ().edge_to (edge)) 
          ? 1.0 : -1.0;
        const double area = geometry ().edge_area (edge); // [cm^2 S]
        const double J_edge       // [g/cm^2 S/h]
          = J_tertiary[edge] + J_secondary[edge] + J_primary[edge];
        J_sum += in_sign * J_edge * area; // [g/h]
        if (in_sign * J_tertiary[edge] < 0.0)
          {
            std::ostringstream tmp;
            tmp << "J_tertiary[" << edge << "] = " << J_tertiary[edge]
                << ", in_sign = " << in_sign << ", J_above = " << J_above;
            msg.bug (tmp.str ());
          }
        if (in_sign * J_secondary[edge] < 0.0)
          {
            std::ostringstream tmp;
            tmp << "J_secondary[" << edge << "] = " << J_secondary[edge]
                << ", in_sign = " << in_sign << ", J_above = " << J_above;
            msg.bug (tmp.str ());
          }
      }
    J_sum /= geometry ().surface_area (); // [g/cm^2 S/h]
    daisy_approximate (-J_above, J_sum);
  }

  // We set a fixed concentration below lower boundary, if specified.
  std::map<size_t, double> C_border;

  const double C_below = chemical.C_below ();
  if (C_below >= 0.0)
    {
      const std::vector<size_t>& edge_below 
        = geometry ().cell_edges (Geometry::cell_below);
      const size_t edge_below_size = edge_below.size ();

      for (size_t i = 0; i < edge_below_size; i++)
        {
          const size_t edge = edge_below[i];
          C_border[edge] = C_below;
        }
    }

  // Tertiary transport.
  tertiary->solute (geometry (), soil_water, J_tertiary, dt, chemical, msg);

  // Fully adsorbed.
  if (chemical.adsorption ().full ())
    {
      static const symbol solid_name ("immobile transport");
      Treelog::Open nest (msg, solid_name);
      if (!iszero (J_above))
        {
          std::ostringstream tmp;
          tmp << "J_above = " << J_above << ", expected 0 for full sorbtion";
          msg.error (tmp.str ());
        }

      // Secondary "transport".
      std::vector<double> J2 (edge_size, 0.0); // Flux delivered by flow.
      std::vector<double> Mn (cell_size); // New content.
//.........这里部分代码省略.........

void 
SummaryBalance::summarize (Treelog& msg) const
{
  TREELOG_MODEL (msg);

  // We write the summary to a string at first.
  std::ostringstream tmp;
  tmp.precision (precision);
  tmp.flags (std::ios::right | std::ios::fixed);
  if (description.name () != default_description)
    tmp << description << "\n\n";

  // Find width of tags.
  const std::string total_title = "Balance (= In - Out - Increase)";
  const std::string content_title = "Total increase in content";
  size_t max_size = std::max (total_title.size (), content_title.size ());
  for (unsigned int i = 0; i < fetch.size (); i++)
    max_size = std::max (max_size, fetch[i]->name_size ());

  // Find width and total values
  int max_digits = 0;
  const double total_input = find_total (input, max_digits);
  const double total_output = find_total (output, max_digits);
  const double total_content = find_total (content, max_digits);
  const double total = total_input - total_output - total_content;
  max_digits = std::max (max_digits, FetchPretty::width (total_input));
  max_digits = std::max (max_digits, FetchPretty::width (total_output));
  max_digits = std::max (max_digits, FetchPretty::width (total_content));
  max_digits = std::max (max_digits, FetchPretty::width (total));

  // Find total width.
  const int width = max_digits + (precision > 0 ? 1 : 0) + precision;

  // Find width of dimensions.
  size_t dim_size = 0;
  for (unsigned int i = 0; i < fetch.size (); i++)
    dim_size = std::max (dim_size, fetch[i]->dimension ().name ().size ());

  // Print all entries.
  symbol shared_dim = Attribute::User ();
  if (input.size () > 0)
    {
      const symbol dim = print_entries (tmp, input, max_size, width);
      print_balance (tmp, "Total input", total_input, dim,
                     dim_size, max_size, width);
      shared_dim = dim;
      tmp << "\n";
    }

  if (output.size () > 0)
    {
      const symbol dim = print_entries (tmp, output, max_size, width);
      print_balance (tmp, "Total output", total_output, dim,
                     dim_size, max_size, width);
      if (shared_dim == Attribute::User ()) 
        shared_dim = dim;
      else if (dim != shared_dim)
        shared_dim = Attribute::Unknown ();
      tmp << "\n";
    }

  if (content.size () > 0)
    {
      const symbol dim = print_entries (tmp, content, max_size, width);
      print_balance (tmp, content_title, total_content, dim,
                     dim_size, max_size, width);
      if (shared_dim == Attribute::User ()) 
        shared_dim = dim;
      else if (dim != shared_dim)
        shared_dim = Attribute::Unknown ();
      tmp << "\n";
    }

  print_balance (tmp, total_title, total, 
                 shared_dim, dim_size, max_size, width);
  tmp << std::string (max_size + 3, ' ') << std::string (width, '=');

  // Where?
  if (file == "")
    msg.message (tmp.str ());
  else
    { 
      std::ofstream out (file.name ().c_str ());
      out << tmp.str ();
      if (! out.good ())
        msg.error ("Could not write to '" + file + "'");
    } 
}

void 
UZRectMollerup::tick (const GeometryRect& geo,
		      const std::vector<size_t>& drain_cell,
		      const double drain_water_level,
		      const Soil& soil, 
		      SoilWater& soil_water, const SoilHeat& soil_heat,
		      const Surface& surface, const Groundwater& groundwater,
		      const double dt, Treelog& msg)

{
  daisy_assert (K_average.get ());
  const size_t edge_size = geo.edge_size (); // number of edges 
  const size_t cell_size = geo.cell_size (); // number of cells 

  // Insert magic here.
  
  ublas::vector<double> Theta (cell_size); // water content 
  ublas::vector<double> Theta_previous (cell_size); // at start of small t-step
  ublas::vector<double> h (cell_size); // matrix pressure
  ublas::vector<double> h_previous (cell_size); // at start of small timestep
  ublas::vector<double> h_ice (cell_size); // 
  ublas::vector<double> S (cell_size); // sink term
  ublas::vector<double> S_vol (cell_size); // sink term
#ifdef TEST_OM_DEN_ER_BRUGT
  ublas::vector<double> S_macro (cell_size);  // sink term
  std::vector<double> S_drain (cell_size, 0.0); // matrix-> macro -> drain flow 
  std::vector<double> S_drain_sum (cell_size, 0.0); // For large timestep
  const std::vector<double> S_matrix (cell_size, 0.0);  // matrix -> macro 
  std::vector<double> S_matrix_sum (cell_size, 0.0); // for large timestep
#endif
  ublas::vector<double> T (cell_size); // temperature 
  ublas::vector<double> Kold (edge_size); // old hydraulic conductivity
  ublas::vector<double> Ksum (edge_size); // Hansen hydraulic conductivity
  ublas::vector<double> Kcell (cell_size); // hydraulic conductivity
  ublas::vector<double> Kold_cell (cell_size); // old hydraulic conductivity
  ublas::vector<double> Ksum_cell (cell_size); // Hansen hydraulic conductivity
  ublas::vector<double> h_lysimeter (cell_size);
  std::vector<bool> active_lysimeter (cell_size);
  const std::vector<size_t>& edge_above = geo.cell_edges (Geometry::cell_above);
  const size_t edge_above_size = edge_above.size ();
  ublas::vector<double> remaining_water (edge_above_size);
  std::vector<bool> drain_cell_on (drain_cell.size (),false); 
  

  for (size_t i = 0; i < edge_above_size; i++)
    {
      const size_t edge = edge_above[i];
      remaining_water (i) = surface.h_top (geo, edge);
    }
  ublas::vector<double> q;	// Accumulated flux
  q = ublas::zero_vector<double> (edge_size);
  ublas::vector<double> dq (edge_size); // Flux in small timestep.
  dq = ublas::zero_vector<double> (edge_size);

  //Make Qmat area diagonal matrix 
  //Note: This only needs to be calculated once... 
  ublas::banded_matrix<double> Qmat (cell_size, cell_size, 0, 0);
  for (int c = 0; c < cell_size; c++)
    Qmat (c, c) = geo.cell_volume (c);
 
  // make vectors 
  for (size_t cell = 0; cell != cell_size ; ++cell) 
    {				
      Theta (cell) = soil_water.Theta (cell);
      h (cell) =  soil_water.h (cell);
      h_ice (cell) = soil_water.h_ice (cell);
      S (cell) =  soil_water.S_sum (cell);
      S_vol (cell) = S (cell) * geo.cell_volume (cell);
      if (use_forced_T)
	T (cell) = forced_T;
      else 
	T (cell) = soil_heat.T (cell); 
      h_lysimeter (cell) = geo.zplus (cell) - geo.cell_z (cell);
    }

  // Remember old value.
  Theta_error = Theta;

  // Start time loop 
  double time_left = dt;	// How much of the large time step left.
  double ddt = dt;		// We start with small == large time step.
  int number_of_time_step_reductions = 0;
  int iterations_with_this_time_step = 0;
  

  int n_small_time_steps = 0;
  
  while (time_left > 0.0)
    {
      if (ddt > time_left)
	ddt = time_left;

      std::ostringstream tmp_ddt;
      tmp_ddt << "Time t = " << (dt - time_left) 
              << "; ddt = " << ddt
              << "; steps " << n_small_time_steps 
              << "; time left = " << time_left;
      Treelog::Open nest (msg, tmp_ddt.str ());

      if (n_small_time_steps > 0
//.........这里部分代码省略.........

 void doIt (Daisy& daisy, const Scope&, Treelog& msg)
 { 
   msg.message ("Sowing " + crop->type_name ());      
   daisy.field ().sow (metalib, *crop, row_width, row_pos, seed, 
                     daisy.time (), msg); 
 }

  double solve_2D (std::vector<Iterative::PointValue>& obs,
                 Treelog& msg)
  {
    // Find best fit.
    struct ToMinimize : Iterative::PointFunction
    {
      const std::vector<Iterative::PointValue>& obs;
      GP2Dfun& fun;
      const double fixed_SoilDepth; // [cm]
      const bool find_SoilDepth;
      const int debug;
      Treelog& msg;

      double value (const Iterative::Point& p) const
      {
        Treelog::Open nest (msg, "minimize");

        if (find_SoilDepth)
          daisy_assert (p.size () == 4);
        else
          daisy_assert (p.size () == 3);
        const double CropDepth = p[0];
        const double CropWidth = p[1];
        const double WRoot = p[2];
        const double SoilDepth = find_SoilDepth ? p[3] : fixed_SoilDepth;

        if (debug > 0)
          {
            std::ostringstream tmp;
            tmp << "CropDepth = " << CropDepth << " cm\n"
                << "CropWidth = " << CropWidth << " cm\n"
                << "WRoot = " << (0.01 * WRoot)  << " Mg DM/ha\n";
            if (find_SoilDepth)
              tmp << "SoilDepth = " << SoilDepth << " cm";
            msg.message (tmp.str ());
          }
       
        // Restrictions.
        const double LARGE_NUMBER = 42.42e42;
        if (CropDepth <= 0)
          return LARGE_NUMBER;
        if (CropWidth <= 0)
          return LARGE_NUMBER;
        if (WRoot <= 0)
          return LARGE_NUMBER;
        if (find_SoilDepth && SoilDepth <= 0)
          return LARGE_NUMBER;

        bool ok = fun.root.set_dynamic (SoilDepth, CropDepth, CropWidth, WRoot,
                                        debug, msg);
        if (!ok)
          return LARGE_NUMBER;

        const double Rsqr = Iterative::RSquared (obs, fun);
        if (debug > 0)
          {
            std::ostringstream tmp;
            tmp << "R^2 = " << Rsqr;
            msg.message (tmp.str ());
          }

        return -Rsqr;
      }

      ToMinimize (const std::vector<Iterative::PointValue>& o, GP2Dfun& f,
                  const double sd,
                  const int d, 
                  Treelog& m)
        : obs (o),
          fun (f),
          fixed_SoilDepth (sd),
          find_SoilDepth (!std::isfinite (fixed_SoilDepth)),
          debug (d),
          msg (m)
      { }
    };
    ToMinimize to_minimize (obs, gp2d, SoilDepth, debug, msg);

    // Initialial guess.
    const double default_CropDepth = 70;  // [cm]
    const double default_CropWidth = 100; // [cm]
    const double default_WRoot = 50; // 150 [g DM/m^2] = 1.5 [Mg DM/ha]
    const double default_SoilDepth = 150; // [cm]
    Iterative::Point start;
    start.push_back (default_CropDepth);
    start.push_back (default_CropWidth);
    start.push_back (default_WRoot);
    if (!std::isfinite (SoilDepth))
      start.push_back (default_SoilDepth);  // [cm]
    daisy_assert (start.size () == 4 || start.size () == 3);
    const double epsilon = 0.01;
    const size_t min_iter = 10000;
    const size_t max_iter = 300000;
    Iterative::Point result;
    const bool solved = Iterative::NelderMead (min_iter, max_iter, epsilon,
                                               to_minimize, start, result,
                                               Treelog::null ());
    const double Rsqr = -to_minimize.value (result);
    
    store ("2D R^2", Rsqr, solved ? "(solved)" : "(no solution)");
//.........这里部分代码省略.........

  double solve_1D (std::vector<Iterative::PointValue>& obs,
                   Treelog& msg)
  {
    // Find best fit.
    struct ToMinimize : Iterative::PointFunction
    {
      const std::vector<Iterative::PointValue>& obs;
      GP1Dfun& fun;
      const double fixed_SoilDepth; // [cm]
      const bool find_SoilDepth;
      const int debug;
      Treelog& msg;

      double value (const Iterative::Point& p) const
      {
        Treelog::Open nest (msg, "minimize");

        if (find_SoilDepth)
          daisy_assert (p.size () == 3);
        else
          daisy_assert (p.size () == 2);
        const double CropDepth = p[0];
        const double WRoot = p[1];
        const double SoilDepth = find_SoilDepth ? p[2] : fixed_SoilDepth;

        if (debug > 0)
          {
            std::ostringstream tmp;
            tmp << "CropDepth = " << CropDepth << " cm\n"
                << "WRoot = " << (0.01 * WRoot)  << " Mg DM/ha\n";
            if (find_SoilDepth)
              tmp << "SoilDepth = " << SoilDepth << " cm";
            msg.message (tmp.str ());
          }
       
        // Restrictions.
        const double LARGE_NUMBER = 42.42e42;
        if (CropDepth <= 0)
          return LARGE_NUMBER;
        if (WRoot <= 0)
          return LARGE_NUMBER;
        if (find_SoilDepth && SoilDepth <= 0)
          return LARGE_NUMBER;

        bool ok = fun.root.set_dynamic (SoilDepth, CropDepth, WRoot,
                                        debug, msg);
        if (!ok)
          return LARGE_NUMBER;

        const double Rsqr = Iterative::RSquared (obs, fun);
        if (debug > 0)
          {
            std::ostringstream tmp;
            tmp << "R^2 = " << Rsqr;
            msg.message (tmp.str ());
          }

        return -Rsqr;
      }

      ToMinimize (const std::vector<Iterative::PointValue>& o, GP1Dfun& f,
                  const double sd,
                  const int d, 
                  Treelog& m)
        : obs (o),
          fun (f),
          fixed_SoilDepth (sd),
          find_SoilDepth (!std::isfinite (fixed_SoilDepth)),
          debug (d),
          msg (m)
      { }
    };
    ToMinimize to_minimize (obs, gp1d, SoilDepth, debug, msg);

    // Initialial guess.
    const double default_CropDepth = 70;  // [cm]
    const double default_WRoot = 50; // 150 [g DM/m^2] = 1.5 [Mg DM/ha]
    const double default_SoilDepth = 150; // [cm]
    Iterative::Point start;
    start.push_back (default_CropDepth);
    start.push_back (default_WRoot);
    if (!std::isfinite (SoilDepth))
      start.push_back (default_SoilDepth);  // [cm]
    daisy_assert (start.size () == 3 || start.size () == 2);
    const double epsilon = 0.01;
    const size_t min_iter = 10000;
    const size_t max_iter = 300000;
    Iterative::Point result;
    const bool solved = Iterative::NelderMead (min_iter, max_iter, epsilon,
                                               to_minimize, start, result,
                                               Treelog::null ());
    const double Rsqr = -to_minimize.value (result);
    
    store ("1D R^2", Rsqr, solved ? "(solved)" : "(no solution)");
    store ("1D CropDepth", result[0], "cm");
    store ("1D WRoot", (0.01 * result[1]) , "Mg DM/ha");;
    if (result.size () == 3)
      store ("1D SoilDepth", result[2], "cm");;
        
    store ("1D L0", gp1d.root.L0, "cm/cm^3");
//.........这里部分代码省略.........

double
PhotoFarquhar::assimilate (const Units& units,
                           const double ABA, const double psi_c,
                           const double ec /* Canopy Vapour Pressure [Pa] */, 
                           const double gbw_ms /* Boundary layer [m/s] */,
			   const double CO2_atm, const double O2_atm, 
                           const double Ptot /* [Pa] */,
			   const double, const double Tc, const double Tl,
                           const double cropN,
			   const std::vector<double>& PAR, 
			   const std::vector<double>& PAR_height,
			   const double PAR_LAI,
			   const std::vector<double>& fraction,
                           const double,
			   CanopyStandard& canopy,
			   Phenology& development,
			   Treelog& msg) 
{
  const double h_x = std::fabs (psi_c) * 1.0e-4  /* [MPa/cm] */; // MPa
  // sugar production [gCH2O/m2/h] by canopy photosynthesis.
  const PLF& LAIvsH = canopy.LAIvsH;
  const double DS = development.DS;

  
  // One crop: daisy_assert (approximate (canopy.CAI, bioclimate.CAI ()));
  if (!approximate (LAIvsH (canopy.Height), canopy.CAI))
    {
      std::ostringstream tmp;
      tmp << "Bug: CAI below top: " << LAIvsH (canopy.Height)
	  << " Total CAI: " << canopy.CAI << "\n";
      canopy.CanopyStructure (DS);
      tmp << "Adjusted: CAI below top: " << LAIvsH (canopy.Height)
	  << " Total CAI: " << canopy.CAI;
      msg.error (tmp.str ());
    }
 
  // CAI (total)  below the current leaf layer.
  double prevLA = LAIvsH (PAR_height[0]); 

  // Assimilate produced by canopy photosynthesis
  double Ass_ = 0.0;

  // Accumulated CAI, for testing purposes.
  double accCAI =0.0;

  // Number of computational intervals in the canopy.
  const int No = PAR.size () - 1;
  daisy_assert (No > 0);
  daisy_assert (No == PAR_height.size () - 1);
  
  // N-distribution and photosynthetical capacity 
  std::vector<double> rubisco_Ndist (No, 0.0);
  std::vector<double> crop_Vm_total (No, 0.0);

  // Photosynthetic capacity (for logging)
  while (Vm_vector.size () < No)
    Vm_vector.push_back (0.0);
  // Potential electron transport rate (for logging)
  while (Jm_vector.size () < No)
    Jm_vector.push_back (0.0);
  // Photosynthetic N-leaf distribution (for logging)
  while (Nleaf_vector.size () < No)
    Nleaf_vector.push_back (0.0);
  // Brutto assimilate production (for logging)
  while (Ass_vector.size () < No)
    Ass_vector.push_back (0.0);
  // LAI (for logging)
  while (LAI_vector.size () < No)
    LAI_vector.push_back (0.0);

  rubiscoNdist->rubiscoN_distribution (units,
                                       PAR_height, prevLA, DS,
                                       rubisco_Ndist/*[mol/m²leaf]*/, 
				       cropN /*[g/m²area]*/, msg);
  crop_Vmax_total (rubisco_Ndist, crop_Vm_total);  


  // Net photosynthesis (for logging)
  while (pn_vector.size () < No)
    pn_vector.push_back (0.0);//

  // Stomata CO2 pressure (for logging)
  while (ci_vector.size () < No)
    ci_vector.push_back (0.0);//[Pa]

  // Leaf surface CO2 pressure  (for logging)
  while (cs_vector.size () < No)
    cs_vector.push_back (0.0);//

  // Leaf surface relative humidity (for logging)
  while (hs_vector.size () < No)
    hs_vector.push_back (0.0);//[]

  // Stomata conductance (for logging)
  while (gs_vector.size () < No)
    gs_vector.push_back (0.0);//[m/s]
     
  // Photosynthetic effect of Xylem ABA and crown water potential.
  Gamma = Arrhenius (Gamma25, Ea_Gamma, Tl); // [Pa]
  const double estar = FAO::SaturationVapourPressure (Tl); // [Pa]
//.........这里部分代码省略.........

bool 
Frame::Implementation::check (const Metalib& metalib, const Frame& frame,
                              const Type& type,
                              const symbol key, Treelog& msg) const
{ 
  // Has value?
  if (!frame.check (key))
    {
      if (type.is_mandatory ())
        {
          msg.error (key + " is missing");
          return false;
        }
      return true;
    }

  // Now, the real checks.
  bool ok = true;

  // Genenic check.
  vcheck_map::const_iterator i = val_checks.find (key);
  if (i != val_checks.end ())
    {
      if (!(*i).second->verify (metalib, frame, key, msg))
        ok = false;;
    }
      
  // Spcial handling of various types.
  switch (type.type ())
    {
    case Attribute::Number:
      // This should already be checked by the file parser, but you
      // never know... Well, theoretically the alist could come from
      // another source (like one of the many other parsers :/), or
      // be one of the build in ones.
      if (type.size () != Attribute::Singleton)
        {
          const std::vector<double>& array = frame.number_sequence (key);
          for (size_t i = 0; i < array.size (); i++)
            {
              std::ostringstream tmp;
              tmp << key << "[" << i << "]";
              Treelog::Open nest (msg, tmp.str ());
              if (!type.verify (array[i], msg))
                ok = false;
            }
        }
      else 
        {
          Treelog::Open nest (msg, key);
          if (!type.verify (frame.number (key), msg))
            ok = false;
        }
      break;

    case Attribute::Model:
      if (type.size () != Attribute::Singleton)
        {
          const ::Library& lib = metalib.library (type.component ());
          const std::vector<boost::shared_ptr<const FrameModel>/**/>& seq 
            = frame.model_sequence (key);
          int j_index = 0;
          for (std::vector<boost::shared_ptr<const FrameModel>/**/>::const_iterator j
                 = seq.begin ();
               j != seq.end ();
               j++)
            {
              std::ostringstream tmp;
              tmp << key << "[" << j_index << "]: ";
              j_index++;
              const FrameModel& al = **j;
              if (!lib.check (al.type_name ()))
                {
                  tmp << "Unknown library member '"
                      << al.type_name () << "'";
                  msg.error (tmp.str ());
                  ok = false;
                }
              else 
                {
                  tmp << al.type_name ();
                  Treelog::Open nest (msg, tmp.str ());
                  if (!al.check (metalib, msg))
                    ok = false;
                }
            }
        }
      else 
        {
          const FrameModel& al = frame.model (key);
          Treelog::Open nest (msg, key + ": " + al.type_name ());
          if (!al.check (metalib, msg))
            ok = false;
        }
      break;

    case Attribute::Submodel:
      if (type.size () != Attribute::Singleton)
        {
          daisy_assert (frame.type_size (key) != Attribute::Singleton);
//.........这里部分代码省略.........

 // Simulation.
 void doIt (Daisy& daisy, const Scope&, Treelog& out)
 {
   out.message ("Adjusting surface detention capacity");
   daisy.field ().set_surface_detention_capacity (height);
 }