C++ Soil类代码示例

void
Ridge::Implementation::update_water (const Geometry1D& geo,
                                     const Soil& soil,
				     const std::vector<double>& S_,
				     std::vector<double>& h_,
				     std::vector<double>& Theta_,
				     std::vector<double>& q,
				     const std::vector<double>& q_p,
                                     const double dt)
{
  const double E = -(q[last_cell + 1] + q_p[last_cell + 1]);
  Theta = (I - E) * dt;
  for (int i = 0; i <= last_cell; i++)
    Theta += (Theta_[i] - S_[i] * dt) * geo.dz (i);
  Theta /= dz;
  const double Theta_sat = soil.Theta (0, 0.0, 0.0);
  daisy_assert (Theta < Theta_sat);
  h = soil.h (0, Theta);

  q[0] = -I;
  for (int i = 0; i <= last_cell; i++)
    {
      q[i+1] = q[i] + geo.dz (i) * (S_[i] + (Theta - Theta_[i]) / dt ) 
	- q_p[i+1];
      Theta_[i] = Theta;
      h_[i] = h;
      if (!approximate (soil.h (i, Theta), h))
	{
	  daisy_assert (i > 0);
	  throw ("Soil hydraulic paramteres change in ridge area");
	}
    }
  daisy_assert (approximate (E, -(q[last_cell + 1] + q_p[last_cell + 1])));
}

double
SoilWater::overflow (const Geometry& geo, 
                     const Soil& soil, const SoilHeat& soil_heat, 
                     Treelog& msg)
{
  const size_t cell_size = geo.cell_size ();
  double extra = 0.0;           // [cm^3]
  for (size_t c = 0; c < cell_size; c++)
    {
      const double h_sat = 0.0;
      const double Theta_sat = soil.Theta (c, h_sat, h_ice (c));
      const double Theta_extra = std::max (Theta_[c] - Theta_sat, 0.0);
      Theta_[c] -= Theta_extra;
      tillage_[c] -= Theta_extra;
      extra += Theta_extra * geo.cell_volume (c);
      const double new_h = soil.h (c, Theta_[c]);
      if (h_[c] < 0.0 || new_h < 0)
        h_[c] = new_h;
    }
  tick_after (geo, soil, soil_heat, false, msg);

  extra /= geo.surface_area (); // [cm]
  extra *= 10.0;                // [mm]
  return extra;
}

std::vector<double> 
Movement1D::default_heat (const Soil& soil, 
                          const Time& time, const Weather& weather)
{
  // Fetch average temperatur.
  const double rad_per_day = 2.0 * M_PI / 365.0;

  // Calculate delay.
  const double pF_2_0 = -100.0;
  double k = 0;
  double C = 0;

  std::vector<double> T;
  const size_t cell_size = geo->cell_size ();

  for (unsigned int i = 0; i < cell_size; i++)
    {
      const double Theta_pF_2_0 = soil.Theta (i, pF_2_0, 0.0);
      k += geo->dz (i) * soil.heat_conductivity (i, Theta_pF_2_0, 0.0);
      C += geo->dz (i) * soil.heat_capacity (i, Theta_pF_2_0, 0.0);
      const double a = k / C;
      delay = geo->zplus (i) / sqrt (24.0 * 2.0 * a / rad_per_day);

      T.push_back (bottom_heat (time, weather));
    }
  daisy_assert (T.size () == cell_size);
  return T;
}

double 
UZRectMollerup::find_K_edge (const Soil& soil, const Geometry& geo, 
                             const size_t e,
                             const ublas::vector<double>& h, 
                             const ublas::vector<double>& h_ice, 
                             const ublas::vector<double>& h_old, 
                             const ublas::vector<double>& T) const
{
  const double anisotropy = soil.anisotropy_edge (e);
  const int from = geo.edge_from (e);
  const int to = geo.edge_to (e);

  // External edges.
  if (!geo.cell_is_internal (from))
    return soil.K (to, h (to), h_ice (to), T (to)) * anisotropy;

  if (!geo.cell_is_internal (to))
    return soil.K (from, h (from), h_ice (from), T (from)) * anisotropy;
  
  // Internal edges.
  const double K_from = soil.K (from, h (from), h_ice (from), T (from));
  const double K_to = soil.K (to, h (to), h_ice (to), T (to));
  return  K_average->average (soil, geo, e, 
                              K_from, h (from), h_ice (from), h_old (from), T (from),
                              K_to, h (to), h_ice (to), h_old (from), T (to)) * anisotropy;
}

double
DrainLateral::K_to_pipes (const unsigned int i, 
                            const Soil& soil, 
                            const SoilHeat& soil_heat) const
{
  if (K_to_pipes_ < 0)
    return soil.K (i, 0.0, 0.0, soil_heat.T (i))
      * soil.anisotropy_cell (i);
  return K_to_pipes_;
}

void 
ReactionNitrification::tick_soil (const Units&, const Geometry& geo,
                                  const Soil& soil, const SoilWater& soil_water,
                                  const SoilHeat& soil_heat,
                                  const OrganicMatter& organic_matter, 
                                  Chemistry& chemistry, 
                                  const double /* dt */, Treelog&)
{
  const size_t cell_size = geo.cell_size ();
  const std::vector<bool> active = organic_matter.active (); 
  Chemical& soil_NO3 = chemistry.find (Chemical::NO3 ());
  Chemical& soil_NH4 = chemistry.find (Chemical::NH4 ());

  daisy_assert (NH4.size () == cell_size);
  daisy_assert (N2O.size () == cell_size);
  daisy_assert (NO3.size () == cell_size);

  for (size_t i = 0; i < cell_size; i++)
    {
      if (active[i])
        soil.nitrification (i, 
                            soil_NH4.M_primary (i), 
			    soil_NH4.C_primary (i), 
                            soil_water.h (i), soil_heat.T (i),
                            NH4[i], N2O[i], NO3[i]);
      else
        NH4[i] = N2O[i] = NO3[i] = 0.0;        
    }


  soil_NH4.add_to_transform_sink (NH4);
  soil_NO3.add_to_transform_source (NO3);
}

double
SoilWater::MaxExfiltration (const Geometry& geo, const size_t edge,
                            const Soil& soil, const double T) const
{
  const size_t n = geo.edge_other (edge, Geometry::cell_above);
  const double h0 = h (n);
  const double K0 = soil.K (n, h0, h_ice (n), T);
  if (max_exfiltration_gradient > 0.0)
    return K0 * max_exfiltration_gradient;
  const double Cw2 = soil.Cw2 (n, h0);
  const double Theta0 = Theta (n);
  const double Theta_surf = soil.Theta_res (n);
  const double delta_Theta = Theta0 - Theta_surf;
  const double z0 = geo.cell_z (n);
  // Darcy formulated for Theta between middle of node and soil surface.
  return - (K0 / Cw2) * (delta_Theta / z0);
}

void
SoilWater::tick_ice (const Geometry& geo, const Soil& soil, 
                     const double dt, Treelog& msg)
{
  TREELOG_SUBMODEL (msg, "SoilWater");

  const size_t cell_size = geo.cell_size ();

  // Ice first.
  for (size_t i = 0; i < cell_size; i++)
    {
      const double porosity = soil.Theta (i, 0.0, 0.0);
      const double Theta_res = soil.Theta_res (i);

      X_ice_[i] -= S_ice_ice[i] * dt;

      // Move extra ice to buffer.
      const double total_ice = X_ice_[i] + X_ice_buffer_[i];
      if (total_ice > 0.0)
        {
          if (Theta_[i] < Theta_res)
            {
              std::ostringstream tmp;
              tmp << "Theta[" << i << "] = " << Theta_[i]
                  << ", less than Theta_res = " << Theta_res;
              daisy_warning (tmp.str ());
            }
          const double Theta_lim = std::max (Theta_res, 
                                             Theta_[i] - S_ice_water_[i] * dt);

          const double available_space = porosity - Theta_lim - 1e-9;
          X_ice_[i] = std::min (available_space, total_ice);
          X_ice_buffer_[i] = total_ice - X_ice_[i];
        }
      else
        {
          X_ice_[i] = 0.0;
          X_ice_buffer_[i] = total_ice;
        }
      daisy_approximate (X_ice_[i] + X_ice_buffer_[i], total_ice);

      // Update ice pressure.
      h_ice_[i] = soil.h (i, porosity - X_ice_[i]);
    }
}

void 
Ridge::Implementation::initialize (const Geometry1D& geo,
                                   const Soil& soil,
                                   const SoilWater& soil_water)
{
  // Find values depending on soil numerics.
  last_cell = geo.interval_plus (lowest);
  daisy_assert (last_cell+1 < soil.size ());
  dz = 0 - geo.zplus (last_cell);
  K_sat_below = soil.K (last_cell+1, 0.0, 0.0, 20.0);

  // Initialize water content.
  Theta = 0.0;
  for (int i = 0; i <= last_cell; i++)
    Theta += soil_water.Theta (i) * geo.dz (i);
  Theta /= dz;
  Theta_pre = Theta;
  daisy_assert (Theta < soil.Theta (0, 0.0, 0.0));
  h = soil.h (0, Theta);
  K = soil.K (0, 0.0, 0.0, 20.0);
}

void
ReactionDenit::initialize (const Units&, const Geometry&,
                           const Soil& soil, const SoilWater&,
                           const SoilHeat&, const Surface&, Treelog&)
{
  const size_t cell_size = soil.size ();

  converted.insert (converted.begin (), cell_size, 0.0);
  converted_fast.insert (converted_fast.begin (), cell_size, 0.0);
  converted_redox.insert (converted_redox.begin (), cell_size, 0.0);
  potential.insert (potential.begin (), cell_size, 0.0);
  potential_fast.insert (potential_fast.begin (), cell_size, 0.0);
}

void
ClayOMBiomod::set_rates (const Soil& soil, const std::vector<SMB*>& smb) const
{ 
  // We always have two SMB pools in BIOMOD.
  daisy_assert (smb.size () == 2);

  // Total death and maintenance.
  const double t_SMB1 = smb[0]->turnover_rate + smb[0]->maintenance;
  const double t_SMB2 = smb[1]->turnover_rate + smb[1]->maintenance;

  // The ratios we want to change.
  const double r_SMB1 = smb[0]->turnover_rate / t_SMB1;
  const double r_SMB2 = smb[1]->turnover_rate / t_SMB2;

  for (unsigned int i  = 0; i < soil.size (); i++)
    {
      // Find modifier.
      const double f = find_f (r_SMB1, r_SMB2, soil.clay (i));
      daisy_assert (f > 0.0);

      // Update turnover rate and maintence.
      const double clay_rate1 = smb[0]->turnover_rate * f;
      daisy_assert (clay_rate1 < 1.0);
      smb[0]->clay_turnover.push_back (clay_rate1);
      daisy_assert (t_SMB1 >= clay_rate1);
      smb[0]->clay_maintenance.push_back (t_SMB1 - clay_rate1);
      const double clay_rate2 = smb[1]->turnover_rate * f;
      daisy_assert (clay_rate2 < 1.0);
      smb[1]->clay_turnover.push_back (clay_rate2);
      daisy_assert (t_SMB2 >= clay_rate2);
      smb[1]->clay_maintenance.push_back (t_SMB2 - clay_rate2);
    }

  for (unsigned int i = 0; i < smb.size (); i++)
    {
      daisy_assert (smb[i]->clay_turnover.size () == soil.size ());
      daisy_assert (smb[i]->clay_maintenance.size () == soil.size ());
    }
}

double 
Movement1D::surface_snow_T (const Soil& soil,
                            const SoilWater& soil_water,
                            const SoilHeat& soil_heat,
                            const double T_snow,
                            const double K_snow,
                            const double dZs) const
{
  // Information about soil.
  const double K_soil 
    = soil.heat_conductivity (0, soil_water.Theta (0),
                              soil_water.X_ice (0)) 
    * 1e-7 * 100.0 / 3600.0; // [erg/cm/h/dg C] -> [W/m/dg C]
  const double Z = -geo->cell_z (0) / 100.0; // [cm] -> [m]
  const double T_soil = soil_heat.T (0); // [dg C]
  daisy_assert (T_soil > -100.0);
  daisy_assert (T_soil < 50.0);

  const double T = (K_soil / Z * T_soil + K_snow / dZs * T_snow) 
    / (K_soil / Z + K_snow / dZs);
  daisy_assert (T > -100.0);
  daisy_assert (T < 50.0);
  return T;
}

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
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
UZlr::tick (Treelog& msg, const GeometryVert& geo,
            const Soil& soil, const SoilHeat& soil_heat,
	    unsigned int first, const Surface& top, 
            const size_t top_edge,
	    unsigned int last, const Groundwater& bottom, 
            const size_t bottom_edge,
	    const std::vector<double>& S,
	    const 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,
            const size_t q_offset,
            std::vector<double>& q_base,
            const double dt)
{
  double *const q = &q_base[q_offset];
  double q_up = 0.0;
  double q_down = 0.0;
  const Surface::top_t top_type = top.top_type (geo, top_edge);

  if (top_type == Surface::soil)
    {
      // We have a forced pressure top, in the form of a ridge system.
      // Since LR only works with flux top, we use Darcy to simulate a
      // flux top between the first cell (with a forced pressure) and
      // the second cell, and then continue calculating with a flux
      // top from the second cell.
      const double dz = geo.cell_z (first) - geo.cell_z (first+1);
      const double dh = (h_old[first] - h_old[first+1]);
      const double K = std::min (soil.K (first, h_old[first], h_ice[first],
                                         soil_heat.T (first)),
                                 soil.K (first, h_old[first+1], h_ice[first+1],
                                         soil_heat.T (first+1)));
      q_up = -K * (dh/dz + 1.0);

      // We can safely ignore S[first], since the ridge system has
      // already incorporated it.
      first++;

      // New upper limit.
      q[first] = q_up;
    }
  else
    {
      // Limit flux by soil capacity.
      const double K_sat = soil.K (first, 0.0, h_ice[first], 
				   soil_heat.T (first));
      daisy_assert (K_sat > 0.0);

      if (top_type == Surface::forced_pressure)
	{
	  const double dz = 0.0 - geo.cell_z (first);
	  const double dh = top.h_top (geo, top_edge) - h_old[first];
	  q_up = q[first] = -K_sat * (dh/dz + 1.0);
	}
      else
        // Limited water or forced flux.
	q_up = q[first] = std::max (top.q_top (geo, top_edge, dt), -K_sat);
    }

  //  Use darcy for upward movement in the top.
  const bool use_darcy = (h_old[first] < h_fc) && (q_up > 0.0);

  // Intermediate cells.
  for (int i = first; i <= last; i++)
    {
      const double z = geo.cell_z (i);
      const double dz = geo.dz (i);
      const double Theta_sat = soil.Theta (i, 0.0, h_ice[i]);
      const double Theta_res = soil.Theta_res (i);
      const double h_min = pF2h (10.0);
      const double Theta_min = soil.Theta (i, h_min, h_ice[i]);
      double Theta_new = Theta_old[i] - q[i] * dt / dz - S[i] * dt;
      if (Theta_new < Theta_min)
        {
          // Extreme dryness.
          q[i+1] = (Theta_min - Theta_new) * dz / dt;
          Theta[i] = Theta_min;
          h[i] = h_min;
          daisy_assert (std::isfinite (h[i]));
          continue;
        }

      daisy_assert (std::isfinite (h_old[i]));
      const double h_new = Theta_new >= Theta_sat 
        ? std::max (h_old[i], 0.0)
        : soil.h (i, Theta_new);
      daisy_assert (std::isfinite (h_new));
      double K_new = soil.K (i, h_new, h_ice[i], soil_heat.T (i));

      // If we have free drainage bottom, we go for field capacity all
      // the way.  Otherwise, we assume groundwater start at the
      // bottom of the last cell, and attempt equilibrium from there.
      // This asumption is correct for lysimeter bottom, adequate for
      // pressure bottom (where groundwater table is in the last
      // cell), and wrong for forced flux (= pipe drained soil) where
      // the groundwater is usually much higher.  Still, it is better
      // than using h_fc.
//.........这里部分代码省略.........