Add constant snow albedo as an option (#235)

* Add constant snow albedo as an option * Add constant snow albedo as an option; use snow density as a conditional dependency; rename snow emissivity * Add the option of asking for constant snow albedo to constructor
amanzi · Feb 14, 2024 · 26799fb · 26799fb
1 parent 5916b34
commit 26799fb
Show file tree

Hide file tree

Showing 4 changed files with 24 additions and 8 deletions.
diff --git a/src/pks/surface_balance/constitutive_relations/land_cover/albedo_threecomponent_evaluator.cc b/src/pks/surface_balance/constitutive_relations/land_cover/albedo_threecomponent_evaluator.cc
@@ -50,10 +50,14 @@ AlbedoThreeComponentEvaluator::AlbedoThreeComponentEvaluator(Teuchos::ParameterL
   // parameters
   a_ice_ = plist_.get<double>("albedo ice [-]", 0.44);
   a_water_ = plist_.get<double>("albedo water [-]", 0.1168);
+  is_constant_snow_albedo_ = plist_.isParameter("albedo snow [-]");
+  if (is_constant_snow_albedo_) {
+    a_snow_ = plist_.get<double>("albedo snow [-]");
+  }
 
   e_ice_ = plist_.get<double>("emissivity ice [-]", 0.98);
   e_water_ = plist_.get<double>("emissivity water [-]", 0.995);
-  e_snow_ = plist_.get<double>("emissivity ground surface [-]", 0.98);
+  e_snow_ = plist_.get<double>("emissivity snow [-]", 0.98);
 }
 
 // Required methods from EvaluatorSecondaryMonotypeCV
@@ -65,7 +69,10 @@ AlbedoThreeComponentEvaluator::Evaluate_(const State& S,
   auto tag = my_keys_.front().second;
 
   // collect dependencies
-  const auto& snow_dens = *S.Get<CompositeVector>(snow_dens_key_, tag).ViewComponent("cell", false);
+  const Epetra_MultiVector* snow_dens = nullptr;
+  if (!is_constant_snow_albedo_) {
+    snow_dens = &(*(S.GetPtr<CompositeVector>(snow_dens_key_, tag)->ViewComponent("cell", false)));
+  }
   const auto& unfrozen_fraction =
     *S.Get<CompositeVector>(unfrozen_fraction_key_, tag).ViewComponent("cell", false);
 
@@ -81,7 +88,7 @@ AlbedoThreeComponentEvaluator::Evaluate_(const State& S,
 
     for (auto c : lc_ids) {
       // albedo of the snow
-      albedo[2][c] = Relations::CalcAlbedoSnow(snow_dens[0][c]);
+      albedo[2][c] = is_constant_snow_albedo_ ? a_snow_ : Relations::CalcAlbedoSnow((*snow_dens)[0][c]); 
 
       // a and e of water
       albedo[1][c] = unfrozen_fraction[0][c] * a_water_ + (1 - unfrozen_fraction[0][c]) * a_ice_;

diff --git a/src/pks/surface_balance/constitutive_relations/land_cover/albedo_threecomponent_evaluator.hh b/src/pks/surface_balance/constitutive_relations/land_cover/albedo_threecomponent_evaluator.hh
@@ -83,7 +83,8 @@ class AlbedoThreeComponentEvaluator : public EvaluatorSecondaryMonotypeCV {
   Key snow_dens_key_;
   Key unfrozen_fraction_key_;
 
-  double a_water_, a_ice_;
+  bool is_constant_snow_albedo_;
+  double a_water_, a_ice_, a_snow_;
   double e_water_, e_ice_, e_snow_;
 
   // this is horrid, because this cannot yet live in state

diff --git a/src/pks/surface_balance/constitutive_relations/land_cover/albedo_twocomponent_evaluator.cc b/src/pks/surface_balance/constitutive_relations/land_cover/albedo_twocomponent_evaluator.cc
@@ -54,10 +54,14 @@ AlbedoTwoComponentEvaluator::AlbedoTwoComponentEvaluator(Teuchos::ParameterList&
   // parameters
   a_ice_ = plist_.get<double>("albedo ice [-]", 0.44);
   a_water_ = plist_.get<double>("albedo water [-]", 0.1168);
+  is_constant_snow_albedo_ = plist_.isParameter("albedo snow [-]");
+  if (is_constant_snow_albedo_) {
+    a_snow_ = plist_.get<double>("albedo snow [-]");
+  }
 
   e_ice_ = plist_.get<double>("emissivity ice [-]", 0.98);
   e_water_ = plist_.get<double>("emissivity water [-]", 0.995);
-  e_snow_ = plist_.get<double>("emissivity ground surface [-]", 0.98);
+  e_snow_ = plist_.get<double>("emissivity snow [-]", 0.98);
 }
 
 // Required methods from EvaluatorSecondaryMonotypeCV
@@ -68,7 +72,10 @@ AlbedoTwoComponentEvaluator::Evaluate_(const State& S, const std::vector<Composi
   auto tag = my_keys_.front().second;
 
   // collect dependencies
-  const auto& snow_dens = *S.Get<CompositeVector>(snow_dens_key_, tag).ViewComponent("cell", false);
+  const Epetra_MultiVector* snow_dens = nullptr;
+  if (!is_constant_snow_albedo_) {
+    snow_dens = &(*(S.GetPtr<CompositeVector>(snow_dens_key_, tag)->ViewComponent("cell", false)));
+  }
   const auto& ponded_depth =
     *S.Get<CompositeVector>(ponded_depth_key_, tag).ViewComponent("cell", false);
   const auto& unfrozen_fraction =
@@ -85,7 +92,7 @@ AlbedoTwoComponentEvaluator::Evaluate_(const State& S, const std::vector<Composi
 
     for (auto c : lc_ids) {
       // albedo of the snow
-      albedo[1][c] = Relations::CalcAlbedoSnow(snow_dens[0][c]);
+      albedo[1][c] = is_constant_snow_albedo_ ? a_snow_ : Relations::CalcAlbedoSnow(*(snow_dens)[0][c]); 
 
       double albedo_water =
         unfrozen_fraction[0][c] * a_water_ + (1 - unfrozen_fraction[0][c]) * a_ice_;

diff --git a/src/pks/surface_balance/constitutive_relations/land_cover/albedo_twocomponent_evaluator.hh b/src/pks/surface_balance/constitutive_relations/land_cover/albedo_twocomponent_evaluator.hh
@@ -83,7 +83,8 @@ class AlbedoTwoComponentEvaluator : public EvaluatorSecondaryMonotypeCV {
   Key albedo_key_, emissivity_key_;
   Key snow_dens_key_, ponded_depth_key_, unfrozen_fraction_key_;
 
-  double a_ice_, a_water_;
+  bool is_constant_snow_albedo_;
+  double a_ice_, a_water_, a_snow_;
   double e_snow_, e_ice_, e_water_;
 
   // this is horrid, because this cannot yet live in state