Adds option to save div q for ELM

Also communicates this info back to ELM for use in global mass
balance.  A few notes:

1. I'm not sure the signs are right here -- is ELM expecting positive
   or negative water to the column/baseflow/runoff?

2. I'm not entirely sure I'm dealing with infiltration correctly.
   There might be a factor of surface area that I'm missing here.

3. I'm not entirely sure that units are right.

Moral of the story, this probably needs some debugging with ELM to
check that this actually closes the water balance correctly.

src/executables/elm_ats_api/elm_ats_driver.cc

@@ -109,6 +109,10 @@ ELM_ATSDriver::ELM_ATSDriver(const Teuchos::RCP<Teuchos::ParameterList>& plist,
   evap_key_ = Keys::readKey(*plist_, domain_surf_, "evaporation", "evaporation");
   trans_key_ = Keys::readKey(*plist_, domain_subsurf_, "transpiration", "transpiration");
 
+  // baseflow and runoff, respectively
+  div_q_key_ = Keys::readKey(*plist_, domain_subsurf_, "divergence of water flux", "divergence_water_flux");
+  div_q_surf_key_ = Keys::readKey(*plist_, domain_surf_, "divergence of water flux", "divergence_water_flux");
+
   // keys for fields used to convert ELM units to ATS units
   surf_mol_dens_key_ = Keys::readKey(*plist_, domain_surf_, "surface molar density", "molar_density_liquid");
   surf_mass_dens_key_ = Keys::readKey(*plist_, domain_surf_, "surface mass density", "mass_density_liquid");
@@ -222,6 +226,11 @@ ELM_ATSDriver::setup()
   requireAtNext(pres_key_, Amanzi::Tags::NEXT, *S_, "flow")
     .SetMesh(mesh_subsurf_)->AddComponent("cell", AmanziMesh::CELL, 1);
 
+  requireAtNext(div_q_key_, Amanzi::Tags::NEXT, *S_)
+    .SetMesh(mesh_subsurf_)->AddComponent("cell", AmanziMesh::CELL, 1);
+  requireAtNext(div_q_surf_key_, Amanzi::Tags::NEXT, *S_)
+    .SetMesh(mesh_surf_)->AddComponent("cell", AmanziMesh::CELL, 1);
+
   Coordinator::setup();
 }
 
@@ -651,7 +660,29 @@ ELM_ATSDriver::get_water_fluxes(double * const surf_subsurf_flx,
       transpiration[j * ncolumns_ + i] = trans[0][cells[j]] * factor;
     }
   }
-  // unclear how to implement net_subsurface_fluxes and net_runon... --ETC
+
+  // net subsurface fluxes, or "baseflow" is given by the integral up the column of the
+  // divergence of ATS's subsurface water flux, minus infiltration
+  if (net_subsurface_fluxes != nullptr) {
+    const auto& divq = S_->Get<CompositeVector>(div_q_key_, Amanzi::Tags::NEXT).ViewComponent("cell", false);
+    for (int i = 0; i != ncolumns_; ++i) {
+      const double mm_per_mol = 1000.0 / surfdens[0][i];
+      const auto& cells = mesh_subsurf_->columns.getCells(i);
+      net_subsurface_fluxes[i] = -infilt[0][i];
+      for (int j = 0; j != ncells_per_col_; ++j) {
+        net_subsurface_fluxes[i] += divq[cells[j]];
+      }
+      net_subsurface_fluxes[i] *= mm_per_mol;
+    }
+  }
+
+  if (net_runon != nullptr) {
+    const auto& divq = S_->Get<CompositeVector>(div_q_surf_key_, Amanzi::Tags::NEXT).ViewComponent("cell", false);
+    for (int i = 0; i != ncolumns_; ++i) {
+      const double mm_per_mol = 1000.0 / surfdens[0][i];
+      net_runon[i] = divq[0][i] * mm_per_mol;
+    }
+  }
 }
 
 

src/pks/flow/overland_pressure.hh

@@ -49,6 +49,12 @@ Solves the diffusion wave equation for overland flow with pressure as a primary
     * `"diffusion`" ``[pde-diffusion-spec]`` The (forward) diffusion operator,
       see PDE_Diffusion_.
 
+    * `"diagnostic: divergence of fluxes`" ``[bool]`` **false** If true, saves
+      the divergence of fluxes as a separate field for diagnostics,
+      visualization, or coupling with other codes (e.g. column-based codes such
+      as ELM, where the vertical integral of this, minus infiltration, becomes
+      baseflow).
+
     * `"diffusion preconditioner`" ``[pde-diffusion-spec]`` **optional** The
       inverse of the diffusion operator.  See PDE_Diffusion_.  Typically this
       is only needed to set Jacobian options, as all others probably should
@@ -301,6 +307,9 @@ class OverlandPressureFlow : public PK_PhysicalBDF_Default {
 
   bool precon_used_;
   bool precon_scaled_;
+  bool diag_div_q_;
+  Key diag_div_q_key_;
+
 
   // boundary condition data
   Teuchos::RCP<Functions::BoundaryFunction> bc_zero_gradient_;

src/pks/flow/overland_pressure_pk.cc

@@ -117,6 +117,11 @@ OverlandPressureFlow::parseParameterList()
   patm_hard_limit_ = plist_->get<bool>("allow no negative ponded depths", false);
   min_vel_ponded_depth_ = plist_->get<double>("min ponded depth for velocity calculation", 1e-2);
   min_tidal_bc_ponded_depth_ = plist_->get<double>("min ponded depth for tidal bc", 0.02);
+
+  // diagnostics
+  diag_div_q_ = plist_->get<bool>("diagnostic: divergence of fluxes", false);
+  if (diag_div_q_)
+    diag_div_q_key_ = Keys::readKey(*plist_, domain_, "divergence of water flux", "divergence_water_flux");
 }
 
 
@@ -344,6 +349,13 @@ OverlandPressureFlow::SetupOverlandFlow_()
     .SetMesh(mesh_)
     ->SetGhosted()
     ->SetComponent("cell", AmanziMesh::Entity_kind::CELL, 3);
+
+  // -- diagnostic for divergence of fluxes
+  if (diag_div_q_) {
+    requireAtNext(diag_div_q_key_, tag_next_, *S_, name_)
+      .SetMesh(mesh_)
+      ->SetComponent("cell", AmanziMesh::Entity_kind::CELL, 1);
+  }
 };
 
 

src/pks/flow/overland_pressure_ti.cc

@@ -63,6 +63,8 @@ OverlandPressureFlow::FunctionalResidual(double t_old,
 
   // diffusion term, treated implicitly
   ApplyDiffusion_(tag_next_, res.ptr());
+  if (diag_div_q_)
+    S_->GetW<CompositeVector>(diag_div_q_key_, tag_next_, name_) = *res;
 
   // more debugging -- write diffusion/flux variables to screen
   vnames = { "z", "h_old", "h_new", "h+z" };

src/pks/flow/richards.hh

@@ -61,6 +61,12 @@ Solves Richards equation:
    * `"diffusion`" ``[pde-diffusion-spec]`` The (forward) diffusion
      operator, see PDE_Diffusion_.
 
+   * `"diagnostic: divergence of fluxes`" ``[bool]`` **false** If true, saves
+     the divergence of fluxes as a separate field for diagnostics,
+     visualization, or coupling with other codes (e.g. column-based codes such
+     as ELM, where the vertical integral of this, minus infiltration, becomes
+     baseflow).
+
    * `"diffusion preconditioner`" ``[pde-diffusion-spec]``
      **optional** The inverse of the diffusion operator.  See
      PDE_Diffusion_.  Typically this is only needed to set Jacobian
@@ -451,8 +457,10 @@ class Richards : public PK_PhysicalBDF_Default {
   Key deform_key_;
   Key depth_key_;
 
-  // debugging control
+  // debugging and diagnostic control
   bool fixed_kr_;
+  bool diag_div_q_;
+  Key diag_div_q_key_;
 
   // -- access methods
   virtual Teuchos::RCP<Operators::Operator>

src/pks/flow/richards_pk.cc

@@ -130,6 +130,11 @@ Richards::parseParameterList()
     ss_primary_key_ = Keys::readKey(*plist_, domain_surf, "pressure", "pressure");
   }
 
+  // diagnostics
+  diag_div_q_ = plist_->get<bool>("diagnostic: divergence of fluxes", false);
+  if (diag_div_q_)
+    diag_div_q_key_ = Keys::readKey(*plist_, domain_, "divergence of water flux", "divergence_water_flux");
+
   // parse inherited lists
   PK_PhysicalBDF_Default::parseParameterList();
 }
@@ -433,6 +438,13 @@ Richards::SetupRichardsFlow_()
     ->SetGhosted()
     ->SetComponent("cell", AmanziMesh::Entity_kind::CELL, 3);
 
+  // -- diagnostic for divergence of fluxes
+  if (diag_div_q_) {
+    requireAtNext(diag_div_q_key_, tag_next_, *S_, name_)
+      .SetMesh(mesh_)
+      ->SetComponent("cell", AmanziMesh::Entity_kind::CELL, 1);
+  }
+
   // Globalization and other timestep control flags
   // -- predictors
   modify_predictor_with_consistent_faces_ =

src/pks/flow/richards_ti.cc

@@ -59,6 +59,9 @@ Richards::FunctionalResidual(double t_old,
 
   // diffusion term, treated implicitly
   ApplyDiffusion_(tag_next_, res.ptr());
+  if (diag_div_q_)
+    S_->GetW<CompositeVector>(diag_div_q_key_, tag_next_, name_) = *res;
+
   // if (vapor_diffusion_) AddVaporDiffusionResidual_(tag_next_, res.ptr());
 
   // more debugging -- write diffusion/flux variables to screen