Merge remote-tracking branch 'origin/master' into jb/elm_api

src/executables/ats_driver.hh

@@ -23,7 +23,7 @@ called `"main`".  That list contains the following required elements:
 
     * `"cycle driver`" ``[coordinator-spec]``  See below.
     * `"mesh`" ``[mesh-typed-spec-list]`` A list of Mesh_ objects.
-    * `"regions`" ``[region-typedsublist-spec-list]`` A list of Region_ objects.
+    * `"regions`" ``[region-typedinline-spec-list]`` A list of Region_ objects.
     * `"visualization`" ``[visualization-spec-list]`` A list of Visualization_ objects.
     * `"observations`" ``[observation-spec-list]`` An list of Observation_ objects.
     * `"checkpoint`" ``[checkpoint-spec]`` A Checkpoint_ spec.

src/executables/ats_registration_files.hh

@@ -22,5 +22,5 @@
 //#include "ats_sediment_transport_registration.hh"
 #include "models_transport_reg.hh"
 #ifdef ALQUIMIA_ENABLED
-#include "pks_chemistry_reg.hh"
+#  include "pks_chemistry_reg.hh"
 #endif

src/executables/coordinator.cc

@@ -207,7 +207,8 @@ Coordinator::initialize()
     S_->set_time(Amanzi::Tags::NEXT, t0_);
 
     for (Amanzi::State::mesh_iterator mesh = S_->mesh_begin(); mesh != S_->mesh_end(); ++mesh) {
-      if (S_->IsDeformableMesh(mesh->first)) { Amanzi::DeformCheckpointMesh(*S_, mesh->first); }
+      if (S_->IsDeformableMesh(mesh->first) && !S_->IsAliasedMesh(mesh->first))
+        Amanzi::DeformCheckpointMesh(*S_, mesh->first);
     }
   }
 
@@ -544,7 +545,8 @@ Coordinator::advance()
 }
 
 
-bool Coordinator::visualize(bool force)
+bool
+Coordinator::visualize(bool force)
 {
   // write visualization if requested
   bool dump = force;
@@ -564,7 +566,8 @@ bool Coordinator::visualize(bool force)
 }
 
 
-bool Coordinator::checkpoint(bool force)
+bool
+Coordinator::checkpoint(bool force)
 {
   int cycle = S_->get_cycle();
   double time = S_->get_time();

src/executables/coordinator.hh

@@ -45,9 +45,9 @@ class Coordinator {
   void report_memory();
 
   bool advance();
-  bool visualize(bool force=false);
-  bool checkpoint(bool force=false);
-  double get_dt(bool after_fail=false);
+  bool visualize(bool force = false);
+  bool checkpoint(bool force = false);
+  double get_dt(bool after_fail = false);
 
  protected:
   void InitializeFromPlist_();

src/pks/energy/energy_bc_factory.hh

@@ -47,7 +47,7 @@ independently specify diffusive fluxes.  This can also be used in cases where
 the mass flux is prescribed to be zero (e.g. bottom boundaries, where this
 might be the geothermal gradient).
 
-Units are in [W m^-2].
+Units are in **[MW m^-2]**, noting the deviation from SI units!
 
 Example:
 
@@ -79,7 +79,7 @@ This boundary condition sets the total flux of energy, from both advection and
 diffusion.  This is not used all that often in real applications, but is common
 for benchmarks or other testing.
 
-Units are in [W m^-2].
+Units are in **[MW m^-2]**, noting the deviation from SI units!
 
 Example:
 

src/pks/flow/constitutive_relations/porosity/compressible_porosity_evaluator.hh

@@ -27,7 +27,7 @@ converge.
 .. _compressible-porosity-evaluator-spec
 .. admonition:: compressible-porosity-evaluator-spec
 
-   * `"compressible porosity model parameters`" ``[compressible-porosity-model-spec-list]``
+   * `"compressible porosity model parameters`" ``[compressible-porosity-standard-model-spec-list]``
 
    KEYS:
 

src/pks/surface_balance/constitutive_relations/land_cover/LandCover.cc

@@ -34,7 +34,6 @@ LandCover::LandCover(Teuchos::ParameterList& plist)
     albedo_ground(plist.get<double>("albedo of bare ground [-]", NAN)),
     emissivity_ground(plist.get<double>("emissivity of bare ground [-]", NAN)),
     albedo_canopy(plist.get<double>("albedo of canopy [-]", NAN)),
-    emissivity_canopy(plist.get<double>("emissivity of canopy [-]", NAN)),
     beers_k_sw(plist.get<double>("Beer's law extinction coefficient, shortwave [-]", NAN)),
     beers_k_lw(plist.get<double>("Beer's law extinction coefficient, longwave [-]", NAN)),
     snow_transition_depth(plist.get<double>("snow transition depth [m]", NAN)),
@@ -121,8 +120,6 @@ checkValid(const std::string& region, const LandCover& lc, const std::string& pa
     throwInvalid(region, "emissivity of bare ground [-]");
   if (parname == "albedo_ground" && std::isnan(lc.albedo_ground))
     throwInvalid(region, "albedo of bare ground [-]");
-  if (parname == "emissivity_canopy" && std::isnan(lc.emissivity_canopy))
-    throwInvalid(region, "emissivity of canopy [-]");
   if (parname == "albedo_canopy" && std::isnan(lc.albedo_canopy))
     throwInvalid(region, "albedo of canopy [-]");
 

src/pks/surface_balance/constitutive_relations/land_cover/LandCover.hh

@@ -75,8 +75,6 @@ same region-based partitioning.
       cover type's bare ground, ranging from [0,1]
     * `"albedo of canopy [-]`" ``[double]``  **NaN** Albedo of the land cover type's
       canopy, ranging from [0,1]
-    * `"emissivity of canopy [-]`" ``[double]``  **NaN** Emissivity of the land
-      cover type's canopy, ranging from [0,1]
 
     * `"Beer's law extinction coefficient, shortwave [-]`"  ``[double]``  **NaN**
     * `"Beer's law extinction coefficient, longwave [-]`" ``[double]`` **NaN**

src/pks/surface_balance/constitutive_relations/land_cover/canopy_radiation_evaluator.cc

@@ -25,23 +25,30 @@ CanopyRadiationEvaluator::CanopyRadiationEvaluator(Teuchos::ParameterList& plist
   my_keys_.clear();
 
   if (Keys::in("shortwave", akey)) {
-    can_down_sw_key_ = Keys::readKey(plist_, domain_canopy_, "canopy downward shortwave radiation", akey);
+    can_down_sw_key_ =
+      Keys::readKey(plist_, domain_canopy_, "canopy downward shortwave radiation", akey);
   } else {
-    can_down_sw_key_ = Keys::readKey(plist_, domain_canopy_, "canopy downward shortwave radiation", "downward_shortwave_radiation");
+    can_down_sw_key_ = Keys::readKey(plist_,
+                                     domain_canopy_,
+                                     "canopy downward shortwave radiation",
+                                     "downward_shortwave_radiation");
   }
   my_keys_.emplace_back(KeyTag{ can_down_sw_key_, tag });
 
   if (Keys::in("longwave", akey)) {
-    can_down_lw_key_ = Keys::readKey(plist_, domain_canopy_, "canopy downward longwave radiation", akey);
+    can_down_lw_key_ =
+      Keys::readKey(plist_, domain_canopy_, "canopy downward longwave radiation", akey);
   } else {
-    can_down_lw_key_ = Keys::readKey(plist_, domain_canopy_, "canopy downward longwave radiation", "downward_longwave_radiation");
+    can_down_lw_key_ = Keys::readKey(
+      plist_, domain_canopy_, "canopy downward longwave radiation", "downward_longwave_radiation");
   }
   my_keys_.emplace_back(KeyTag{ can_down_lw_key_, tag });
 
   if (Keys::in("balance", akey)) {
     rad_bal_can_key_ = Keys::readKey(plist_, domain_canopy_, "canopy radiation balance", akey);
   } else {
-    rad_bal_can_key_ = Keys::readKey(plist_, domain_canopy_, "canopy radiation balance", "radiation_balance");
+    rad_bal_can_key_ =
+      Keys::readKey(plist_, domain_canopy_, "canopy radiation balance", "radiation_balance");
   }
   my_keys_.emplace_back(KeyTag{ rad_bal_can_key_, tag });
 
@@ -85,7 +92,7 @@ CanopyRadiationEvaluator::Evaluate_(const State& S, const std::vector<CompositeV
 {
   Tag tag = my_keys_.front().second;
   Epetra_MultiVector& down_sw = *results[0]->ViewComponent("cell", false);
-  Epetra_MultiVector& down_lw = *results[0]->ViewComponent("cell", false);
+  Epetra_MultiVector& down_lw = *results[1]->ViewComponent("cell", false);
   Epetra_MultiVector& rad_bal_can = *results[2]->ViewComponent("cell", false);
 
   const Epetra_MultiVector& sw_in =
@@ -105,37 +112,38 @@ CanopyRadiationEvaluator::Evaluate_(const State& S, const std::vector<CompositeV
       lc.first, AmanziMesh::Entity_kind::CELL, AmanziMesh::Parallel_type::OWNED, &lc_ids);
 
     for (auto c : lc_ids) {
-      // Beer's law to find attenuation of radiation to surface in sw
-      double sw_atm_surf = Relations::BeersLaw(sw_in[0][c], lc.second.beers_k_sw, lai[0][c]);
-      double sw_atm_can = sw_in[0][c] - sw_atm_surf;
+      // NOTE: emissivity = absorptivity, we use e to notate both
+      // Beer's law to find absorptivity of canopy
+      double e_can_sw = Relations::BeersLawAbsorptivity(lc.second.beers_k_sw, lai[0][c]);
+      double e_can_lw = Relations::BeersLawAbsorptivity(lc.second.beers_k_lw, lai[0][c]);
 
-      // Beer's law to find attenuation of radiation to surface in lw -- note
-      // this should be almost 0 for any LAI
-      double lw_atm_surf = Relations::BeersLaw(lw_in[0][c], lc.second.beers_k_lw, lai[0][c]);
-      double lw_atm_can = lw_in[0][c] - lw_atm_surf;
+      // sw atm to canopy and surface
+      double sw_atm_can = e_can_sw * sw_in[0][c];
+      double sw_atm_surf = sw_in[0][c] - sw_atm_can;
 
-      // smooth between lai = 0 (no canopy = no outgoing longwave) to lai = 1
-      // (lai of 1 approximately indicates the entire grid cell is covered in
-      // leaf area?)
-      double lai_factor = lai[0][c] < 1. ? lai[0][c] : 1.;
+      // lw atm to canopy and surface
+      double lw_atm_can = e_can_lw * lw_in[0][c];
+      double lw_atm_surf = lw_in[0][c] - lw_atm_can;
 
       // black-body radiation for LW out
-      double lw_can = lai_factor * Relations::OutgoingLongwaveRadiation(
-                                     temp_canopy[0][c], lc.second.emissivity_canopy);
-
+      double lw_can = Relations::OutgoingLongwaveRadiation(temp_canopy[0][c], e_can_lw);
       down_sw[0][c] = sw_atm_surf;
       down_lw[0][c] = lw_atm_surf + lw_can;
-      // factor of 2 is for upward and downward longwave emission -- is this right?
+
+      // factor of 2 is for upward and downward longwave emission NOTE: this is
+      // missing surface-to-canopy LW radiation, which is included externally
+      // with another evaluator as it is calculated later in the surface
+      // balance evaluator.
       rad_bal_can[0][c] = (1 - lc.second.albedo_canopy) * sw_atm_can + lw_atm_can - 2 * lw_can;
     }
   }
 }
 
 void
 CanopyRadiationEvaluator::EvaluatePartialDerivative_(const State& S,
-                                                      const Key& wrt_key,
-                                                      const Tag& wrt_tag,
-                                                      const std::vector<CompositeVector*>& results)
+                                                     const Key& wrt_key,
+                                                     const Tag& wrt_tag,
+                                                     const std::vector<CompositeVector*>& results)
 {
   for (const auto& res : results) res->PutScalar(0.);
 }

src/pks/surface_balance/constitutive_relations/land_cover/pet_priestley_taylor_evaluator.cc

@@ -19,7 +19,7 @@ namespace Relations {
 namespace PriestleyTaylor {
 
 // Convert all units to SI
-// input temperature [K], output latent heat [J/kg] 
+// input temperature [K], output latent heat [J/kg]
 double
 latentHeatVaporization_water(double temp_air)
 {
@@ -38,23 +38,23 @@ latentHeatVaporization_snow(double temp_air)
 double
 psychrometricConstant(double lh_vap, double elev)
 {
-  double elev_ft = elev * 3.281;               // ft
-  double lh_vap_calg = lh_vap / 1000 / 4.184;  // cal/g
+  double elev_ft = elev * 3.281;              // ft
+  double lh_vap_calg = lh_vap / 1000 / 4.184; // cal/g
   double psy_kpaC = 1.6286 * (101.3 - (0.003215 * elev_ft)) / lh_vap_calg;
-                                               // Kpa/C
-  return psy_kpaC * 1000;                      // Pa/C
+  // Kpa/C
+  return psy_kpaC * 1000; // Pa/C
 }
 
 // input temperature [K], output slope [Pa/C]
 double
 vaporPressureSlope(double temp_air)
 {
   double tempC = temp_air - 273.15;                            // C
-  double vp_sat = Relations::SaturatedVaporPressure(temp_air); // Pa 
+  double vp_sat = Relations::SaturatedVaporPressure(temp_air); // Pa
   return 4098 * vp_sat / std::pow(tempC + 237.3, 2);           // Pa/C
 }
 
-// input temperature [K], output heat flux [W/m^2] 
+// input temperature [K], output heat flux [W/m^2]
 double
 groundHeatFlux(double temp_ground, double temp_air)
 {

src/pks/surface_balance/constitutive_relations/land_cover/radiation_balance_evaluator.cc

@@ -82,9 +82,8 @@ void
 RadiationBalanceEvaluator::EnsureCompatibility_ToDeps_(State& S)
 {
   if (!compatible_) {
-    land_cover_ =
-      getLandCover(S.ICList().sublist("land cover types"),
-                   { "beers_k_lw", "beers_k_sw", "emissivity_canopy", "albedo_canopy" });
+    land_cover_ = getLandCover(S.ICList().sublist("land cover types"),
+                               { "beers_k_lw", "beers_k_sw", "albedo_canopy" });
 
     for (const auto& dep : dependencies_) {
       // dependencies on same mesh, but some have two
@@ -141,31 +140,37 @@ RadiationBalanceEvaluator::Evaluate_(const State& S, const std::vector<Composite
       lc.first, AmanziMesh::Entity_kind::CELL, AmanziMesh::Parallel_type::OWNED, &lc_ids);
 
     for (auto c : lc_ids) {
-      // Beer's law to find attenuation of radiation to surface in sw
-      double sw_atm_surf = Relations::BeersLaw(sw_in[0][c], lc.second.beers_k_sw, lai[0][c]);
-      double sw_atm_can = sw_in[0][c] - sw_atm_surf;
+      // NOTE: emissivity = absorptivity, we use e to notate both
+      // Beer's law to find absorptivity of canopy
+      double e_can_sw = Relations::BeersLawAbsorptivity(lc.second.beers_k_sw, lai[0][c]);
+      double e_can_lw = Relations::BeersLawAbsorptivity(lc.second.beers_k_lw, lai[0][c]);
 
-      // Beer's law to find attenuation of radiation to surface in lw -- note
-      // this should be almost 0 for any LAI
-      double lw_atm_surf = Relations::BeersLaw(lw_in[0][c], lc.second.beers_k_lw, lai[0][c]);
-      double lw_atm_can = lw_in[0][c] - lw_atm_surf;
+      // sw atm to canopy and surface
+      double sw_atm_can = e_can_sw * sw_in[0][c];
+      double sw_atm_surf = sw_in[0][c] - sw_atm_can;
 
-      // smooth between lai = 0 (no canopy = no outgoing longwave) to lai = 1
-      // (lai of 1 approximately indicates the entire grid cell is covered in
-      // leaf area?)
-      double lai_factor = lai[0][c] < 1. ? lai[0][c] : 1.;
+      // lw atm to canopy and surface
+      double lw_atm_can = e_can_lw * lw_in[0][c];
+      double lw_atm_surf = lw_in[0][c] - lw_atm_can;
 
-      // black-body radiation for LW out
+      // lw out of each layer
       double lw_surf = Relations::OutgoingLongwaveRadiation(temp_surf[0][c], emiss[0][c]);
       double lw_snow = Relations::OutgoingLongwaveRadiation(temp_snow[0][c], emiss[1][c]);
-      double lw_can = lai_factor * Relations::OutgoingLongwaveRadiation(
-                                     temp_canopy[0][c], lc.second.emissivity_canopy);
-
-      rad_bal_surf[0][c] = (1 - albedo[0][c]) * sw_atm_surf + lw_atm_surf + lw_can - lw_surf;
-      rad_bal_snow[0][c] = (1 - albedo[1][c]) * sw_atm_surf + lw_atm_surf + lw_can - lw_snow;
-      rad_bal_can[0][c] = (1 - lc.second.albedo_canopy) * sw_atm_can + lw_atm_can +
-                          lai_factor * area_frac[1][c] * lw_snow +
-                          lai_factor * area_frac[0][c] * lw_surf - 2 * lw_can;
+      double lw_can = Relations::OutgoingLongwaveRadiation(temp_canopy[0][c], e_can_lw);
+
+      // surface connections
+      double lw_down = lw_atm_surf + lw_can;
+      double lw_up_surf = (1 - emiss[0][c]) * lw_down + lw_surf;
+      double lw_up_snow = (1 - emiss[1][c]) * lw_down + lw_snow;
+
+      // radiation balances -- see Figure 4.1 in CLM Tech Note
+      rad_bal_surf[0][c] = (1 - albedo[0][c]) * sw_atm_surf + lw_down - lw_up_surf;
+      rad_bal_snow[0][c] = (1 - albedo[1][c]) * sw_atm_surf + lw_down - lw_up_snow;
+
+      rad_bal_can[0][c] =  (1 - lc.second.albedo_canopy) * sw_atm_can
+        + lw_atm_can - 2 * lw_can
+        + area_frac[0][c] * e_can_lw * lw_up_surf
+        + area_frac[1][c] * e_can_lw * lw_up_snow;
     }
   }
 }

src/pks/surface_balance/constitutive_relations/land_cover/radiation_balance_evaluator.hh

@@ -10,14 +10,19 @@
 //! Evaluates a net radiation balance for ground and canopy.
 /*!
 
-
-
 Here the net radiation is positive for energy inputs to the layer.  Note that
 ground is based on the two-channel (land + snow) while canopy is assumed to be
 a simple, single layer.
 
-Requires the use of LandCover types, for albedo and emissivity of the canopy
-itself, along with Beer's law coefficients.
+Requires the use of LandCover types, for albedo and Beer's law coefficients.
+
+This is combination of CLM v4.5 Tech Note and Beer's law for attenuation of
+radiation absorption.  In particular, long-wave is exactly as Figure 4.1c in CLM
+4.5 Tech Note.  The main difference comes in how absorptivity (which is equal
+to emissivity, epsilon in that document) is defined.  Here we use Beer's law
+which is an exponential decay with LAI.
+
+Unlike CLM 4.5, here we do not split shortwave into direct and diffuse light.
 
 Computes:
 

src/pks/surface_balance/constitutive_relations/land_cover/seb_physics_funcs.cc

@@ -84,12 +84,11 @@ OutgoingLongwaveRadiation(double temp, double emissivity)
 }
 
 double
-BeersLaw(double sw_in, double k_extinction, double lai)
+BeersLawAbsorptivity(double k_extinction, double lai)
 {
-  return sw_in * std::exp(-k_extinction * lai);
+  return 1 - std::exp(-k_extinction * lai);
 }
 
-
 double
 WindFactor(double Us, double Z_Us, double Z_rough, double KB)
 {

src/pks/surface_balance/constitutive_relations/land_cover/seb_physics_funcs.hh

@@ -79,7 +79,7 @@ OutgoingLongwaveRadiation(double temp, double emissivity);
 // Beer's law for radiation attenuation through a single-layer canopy
 // ------------------------------------------------------------------------------------------
 double
-BeersLaw(double sw_in, double k_extinction, double lai);
+BeersLawAbsorptivity(double k_extinction, double lai);
 
 //
 // Wind speed term D_he