changes the algorithm to set a max pc at the plant collar rather than…

… using CLM beta function

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

@@ -51,7 +51,7 @@ SoilPlantFluxFunctor::SoilPlantFluxFunctor(AmanziMesh::Entity_ID sc_,
 double
 SoilPlantFluxFunctor::operator()(double plant_pc) const
 {
-  double beta = computeTranspirationReductionFunction(plant_pc);
+  double beta = 1.;//computeTranspirationReductionFunction(plant_pc);
   double total_trans = 0.;
 
   // root pc is given by plant pc (assumed at the surface) - hydrostatic
@@ -98,9 +98,10 @@ SoilPlantFluxFunctor::computeSoilPlantFlux(double root_pc, AmanziMesh::Entity_ID
 
 
 // all right hand sides
-void
+double
 SoilPlantFluxFunctor::computeSoilPlantFluxes(double plant_pc, Epetra_MultiVector& res) const
 {
+  double total_trans = 0.;
   double root_pc = plant_pc;
   double sa_col = sa[0][sc] * 2.0;
   for (auto c : cells_of_col) {
@@ -110,11 +111,14 @@ SoilPlantFluxFunctor::computeSoilPlantFluxes(double plant_pc, Epetra_MultiVector
     root_pc += rho_g * dz_on_2;
 
     // compute flux
-    res[0][c] = computeSoilPlantFlux(root_pc, c) * sa[0][sc] / cv[0][c];
+    double ltrans = computeSoilPlantFlux(root_pc, c);
+    res[0][c] = ltrans * sa[0][sc] / cv[0][c];
+    total_trans += ltrans;
 
     // bottom half-cell
     root_pc += rho_g * dz_on_2;
   }
+  return total_trans;
 }
 
 
@@ -240,46 +244,53 @@ TranspirationDistributionRelPermEvaluator::Evaluate_(const State& S,
                                     rho_,
                                     g);
 
-          // bracket the root -- linear to the left of 1, log to the right of 1
-          std::pair<double, double> ab;
-          if (func(0.) > 0.) {
-            ab.first = 0.;
-            ab.second = 1.e4;
-            while (func(ab.second) > 0) {
-              ab.first = ab.second;
-              ab.second *= 10;
-              if (ab.second > 1.e14) {
-                Errors::Message msg;
-                msg << "TranspirationDistributionRelPermEvaluator:: Failing to bracket root: "
-                    << "start = " << 0. << " f_start = " << func(0.) << "  stop = " << ab.second
-                    << " f_stop = " << func(ab.second);
-                Exceptions::amanzi_throw(msg);
-              }
-            }
+          // compute the flux at max pc
+          double pc_plant_max = region_lc.second.stomata_closed_capillary_pressure;
+          double total_trans = func.computeSoilPlantFluxes(pc_plant_max, trans_v);
+          if (total_trans <= potential_trans[0][sc]) {
+            // insufficient water
+            plant_pc_v[0][sc] = pc_plant_max;
           } else {
-            ab.second = 0.;
-            ab.first = -1.e4;
-            while (func(ab.first) < 0) {
-              ab.second = ab.first;
-              ab.first *= 10;
-              if (ab.first < -1.e14) {
-                Errors::Message msg;
-                msg << "TranspirationDistributionRelPermEvaluator:: Failing to bracket root: "
-                    << "start = " << 0. << " f_start = " << func(0.) << "  stop = " << ab.first
-                    << " f_stop = " << func(ab.first);
-                Exceptions::amanzi_throw(msg);
+            // bracket the root -- linear to the left of 1, log to the right of 1
+            std::pair<double, double> ab;
+            if (func(0.) > 0.) {
+              ab.first = 0.;
+              ab.second = 1.e4;
+              while (func(ab.second) > 0) {
+                ab.first = ab.second;
+                ab.second *= 10;
+                if (ab.second > 1.e14) {
+                  Errors::Message msg;
+                  msg << "TranspirationDistributionRelPermEvaluator:: Failing to bracket root: "
+                      << "start = " << 0. << " f_start = " << func(0.) << "  stop = " << ab.second
+                      << " f_stop = " << func(ab.second);
+                  Exceptions::amanzi_throw(msg);
+                }
+              }
+            } else {
+              ab.second = 0.;
+              ab.first = -1.e4;
+              while (func(ab.first) < 0) {
+                ab.second = ab.first;
+                ab.first *= 10;
+                if (ab.first < -1.e14) {
+                  Errors::Message msg;
+                  msg << "TranspirationDistributionRelPermEvaluator:: Failing to bracket root: "
+                      << "start = " << 0. << " f_start = " << func(0.) << "  stop = " << ab.first
+                      << " f_stop = " << func(ab.first);
+                  Exceptions::amanzi_throw(msg);
+                }
               }
             }
-          }
 
-          // compute the plant capillary pressure using a root-finder
-          int itrs = nits_;
-          plant_pc_v[0][sc] = Amanzi::Utils::findRootBrent(func, ab.first, ab.second, tol_, &itrs);
-          AMANZI_ASSERT(itrs > 0 && itrs <= nits_);
-
-          // compute the distributed transpiration fluxes for each grid cell
-          func.computeSoilPlantFluxes(plant_pc_v[0][sc], trans_v);
+            // compute the plant capillary pressure using a root-finder
+            int itrs = nits_;
+            plant_pc_v[0][sc] = Amanzi::Utils::findRootBrent(func, ab.first, ab.second, tol_, &itrs);
+            AMANZI_ASSERT(itrs > 0 && itrs <= nits_);
 
+            // compute the distributed transpiration fluxes for each grid cell
+            func.computeSoilPlantFluxes(plant_pc_v[0][sc], trans_v);
+          }
         } else {
           for (auto c : subsurf_mesh.columns.getCells(sc)) { trans_v[0][c] = 0.; }
         }

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

@@ -124,7 +124,7 @@ struct SoilPlantFluxFunctor {
 
   // right hand side
   double computeSoilPlantFlux(double root_pc, AmanziMesh::Entity_ID c) const;
-  void computeSoilPlantFluxes(double root_pc, Epetra_MultiVector& trans) const;
+  double computeSoilPlantFluxes(double root_pc, Epetra_MultiVector& trans) const;
 
   // beta in left hand side
   double computeTranspirationReductionFunction(double plant_pc) const;