C++ Soil::h方法代码示例

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;
}

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
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 
Ridge::Implementation::tick (const Geometry1D& geo,
                             const Soil& soil, const SoilWater& soil_water,
			     const double external_ponding, const double dt)
{
  // First, we need to find the internal ponding height. 
  // The external ponding assumes a flat surface, we need to find the
  // point (x_pond, z_pond) in the ridge geometry where air, soil and
  // surface water all meet.
  if (external_ponding < 1.0e-5)
    {
      // No ponding.
      x_pond = 0.0;
      z_pond = lowest;
    }
  else if (external_ponding > highest - 1.0e-5)
    {
      // Above top of ridge (ridge geometry doesn't matter then). 
      x_pond = 1.0;
      z_pond = external_ponding;
    }
  else
    {
      // Somewhere in between.  
      //
      // For a given point (x, z), x * z is the amount of space
      // available had all the soil from the ridge been removed.  If
      // you remove the space occupied by the ridge soil, the
      // remaining space is available to free water and air.  So we
      // need to find the point where x * z minus the ridge soil is
      // equal to the amount of water in the pond.  
      //
      // We can't solve this analytically in general, but we can solve
      // it when the ridge geometry is a straight line.  Since the
      // geometry is represented by a PLF, we can find the relevant
      // piece and solve it there.

      // Note that x * z can be negative because our zero point is at
      // the original soil surface, not the bottom of the ridge, this
      // doesn't matter since we are only interested in differences,
      // not absolute number.  

      // Ridge soil until now.
      double integral = 0.0;
      // Last point in the PLF.
      double x0 = 0.0;
      double z0 = lowest;
      
      for (unsigned int i = 1; ; i++)
	{
	  // We already know the answer must lie between two points.
	  daisy_assert (i < z.size ());
	  // New point in the PLF.
	  const double x1 = z.x (i);
	  const double z1 = z.y (i);
	  // Difference.
	  const double delta_z = z1 - z0;
	  const double delta_x = x1 - x0;
	  // Average z height.
	  const double average_z = 0.5 * (z0 + z1);
	  // Ridge soil for this step.
	  const double this_step = delta_x * average_z;
	  // Total space for this step.
	  const double total = x1 * z1;

	  // Check if this interval can contain all the ponded water.
	  if (total - (integral + this_step) >= external_ponding)
	    {
	      // Find slant in this interval.
	      const double slant = delta_z / delta_x;

	      // We now need to find the point (x_p, z_p) where 
	      //     total - (integral + this_step) = external_ponding)
	      // where
	      //     total = x_p * z_p
	      //     this_step = 0.5 * (x_p - x0) * (z0 + z_p)
	      //
	      // Substitute 
	      //     z_p = z0 + slant * (x_p - x0)
	      // and we get
	      //     x_p * (z0 + slant * (x_p - x0))
	      //     - (integral 
	      //        + 0.5 * (x_p - x0) * (z0 + z0 + slant * (x_p - x0))
	      //     = external_ponding
	      // <=>
	      //     x_p * z0 + slant * x_p^2 - slant * x0 * x_p
	      //     - integral
	      //     - (x_p - x0) * ( z0 + 0.5 * slant * (x_p - x0))
	      //     - external_ponding
	      //     = 0
	      // <=>
	      //     x_p * z0 + slant * x_p^2 - slant * x0 * x_p
	      //     - integral 
	      //     - x_p * z0
	      //     - 0.5 * x_p * slant * (x_p - x0)
	      //     + x0 * z0
	      //     + 0.5 * x0 * slant * (x_p - x0)
	      //     - external_ponding
	      //     = 0
	      // <=> 
//.........这里部分代码省略.........

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.
//.........这里部分代码省略.........