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

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