Fixed #38: cell.surface_water_coverage is using cell.surfacewater.get…

…_coverage()

cmf/cmf_core.py

@@ -7333,7 +7333,7 @@ def get_surfacewater(self, *args, **kwargs):
 
     def surfacewater_as_storage(self, *args, **kwargs):
         """
-        surfacewater_as_storage(Cell self)
+        surfacewater_as_storage(Cell self) -> cmf::upslope::surfacewater_ptr
 
         void surfacewater_as_storage()
 
@@ -7422,7 +7422,6 @@ def albedo(self, *args, **kwargs):
         """
         return _cmf_core.Cell_albedo(self, *args, **kwargs)
 
-    surface_amplitude = _swig_property(_cmf_core.Cell_surface_amplitude_get, _cmf_core.Cell_surface_amplitude_set)
 
     def surface_water_coverage(self, *args, **kwargs):
         """
@@ -10585,6 +10584,11 @@ def get_height_function(self, *args, **kwargs):
         return _cmf_core.SurfaceWater_get_height_function(self, *args, **kwargs)
 
 
+    def get_coverage(self, *args, **kwargs):
+        """get_coverage(SurfaceWater self) -> double"""
+        return _cmf_core.SurfaceWater_get_coverage(self, *args, **kwargs)
+
+
     def get_cell(self, *args, **kwargs):
         """
         get_cell(SurfaceWater self) -> Cell
@@ -10610,6 +10614,7 @@ def __repr__(self):
 
     __swig_destroy__ = _cmf_core.delete_SurfaceWater
 SurfaceWater.get_height_function = new_instancemethod(_cmf_core.SurfaceWater_get_height_function, None, SurfaceWater)
+SurfaceWater.get_coverage = new_instancemethod(_cmf_core.SurfaceWater_get_coverage, None, SurfaceWater)
 SurfaceWater.get_cell = new_instancemethod(_cmf_core.SurfaceWater_get_cell, None, SurfaceWater)
 _cmf_core.SurfaceWater_swigregister(SurfaceWater)
 # SurfaceWater end

cmf/cmf_core_src/cmf_wrap.cpp

@@ -45198,6 +45198,7 @@ SWIGINTERN PyObject *_wrap_Cell_surfacewater_as_storage(PyObject *SWIGUNUSEDPARM
   void *argp1 = 0 ;
   int res1 = 0 ;
   PyObject *swig_obj[1] ;
+  cmf::upslope::surfacewater_ptr result;
 
   if (!args) SWIG_fail;
   swig_obj[0] = args;
@@ -45208,15 +45209,15 @@ SWIGINTERN PyObject *_wrap_Cell_surfacewater_as_storage(PyObject *SWIGUNUSEDPARM
   arg1 = reinterpret_cast< cmf::upslope::Cell * >(argp1);
   {
     try {
-      (arg1)->surfacewater_as_storage();
+      result = (arg1)->surfacewater_as_storage();
     } catch (const std::out_of_range& e) {
       SWIG_exception(SWIG_IndexError, e.what());    
     } catch (const std::exception& e) {
       SWIG_exception(SWIG_RuntimeError, e.what());
     }
 
   }
-  resultobj = SWIG_Py_Void();
+  resultobj = SWIG_NewPointerObj((new cmf::upslope::surfacewater_ptr(static_cast< const cmf::upslope::surfacewater_ptr& >(result))), SWIGTYPE_p_std__shared_ptrT_cmf__upslope__SurfaceWater_t, SWIG_POINTER_OWN |  0 );
   return resultobj;
 fail:
   return NULL;
@@ -45723,58 +45724,6 @@ SWIGINTERN PyObject *_wrap_Cell_albedo(PyObject *SWIGUNUSEDPARM(self), PyObject
 }
 
 
-SWIGINTERN PyObject *_wrap_Cell_surface_amplitude_set(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
-  PyObject *resultobj = 0;
-  cmf::upslope::Cell *arg1 = (cmf::upslope::Cell *) 0 ;
-  real arg2 ;
-  void *argp1 = 0 ;
-  int res1 = 0 ;
-  double val2 ;
-  int ecode2 = 0 ;
-  PyObject *swig_obj[2] ;
-
-  if (!SWIG_Python_UnpackTuple(args,"Cell_surface_amplitude_set",2,2,swig_obj)) SWIG_fail;
-  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_cmf__upslope__Cell, 0 |  0 );
-  if (!SWIG_IsOK(res1)) {
-    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "Cell_surface_amplitude_set" "', argument " "1"" of type '" "cmf::upslope::Cell *""'"); 
-  }
-  arg1 = reinterpret_cast< cmf::upslope::Cell * >(argp1);
-  ecode2 = SWIG_AsVal_double(swig_obj[1], &val2);
-  if (!SWIG_IsOK(ecode2)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "Cell_surface_amplitude_set" "', argument " "2"" of type '" "real""'");
-  } 
-  arg2 = static_cast< real >(val2);
-  if (arg1) (arg1)->surface_amplitude = arg2;
-  resultobj = SWIG_Py_Void();
-  return resultobj;
-fail:
-  return NULL;
-}
-
-
-SWIGINTERN PyObject *_wrap_Cell_surface_amplitude_get(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
-  PyObject *resultobj = 0;
-  cmf::upslope::Cell *arg1 = (cmf::upslope::Cell *) 0 ;
-  void *argp1 = 0 ;
-  int res1 = 0 ;
-  PyObject *swig_obj[1] ;
-  real result;
-
-  if (!args) SWIG_fail;
-  swig_obj[0] = args;
-  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_cmf__upslope__Cell, 0 |  0 );
-  if (!SWIG_IsOK(res1)) {
-    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "Cell_surface_amplitude_get" "', argument " "1"" of type '" "cmf::upslope::Cell *""'"); 
-  }
-  arg1 = reinterpret_cast< cmf::upslope::Cell * >(argp1);
-  result = (real) ((arg1)->surface_amplitude);
-  resultobj = SWIG_From_double(static_cast< double >(result));
-  return resultobj;
-fail:
-  return NULL;
-}
-
-
 SWIGINTERN PyObject *_wrap_Cell_surface_water_coverage(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
   PyObject *resultobj = 0;
   cmf::upslope::Cell *arg1 = (cmf::upslope::Cell *) 0 ;
@@ -63029,6 +62978,50 @@ SWIGINTERN PyObject *_wrap_SurfaceWater_get_height_function(PyObject *SWIGUNUSED
 }
 
 
+SWIGINTERN PyObject *_wrap_SurfaceWater_get_coverage(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
+  PyObject *resultobj = 0;
+  cmf::upslope::SurfaceWater *arg1 = (cmf::upslope::SurfaceWater *) 0 ;
+  void *argp1 = 0 ;
+  int res1 = 0 ;
+  std::shared_ptr< cmf::upslope::SurfaceWater const > tempshared1 ;
+  std::shared_ptr< cmf::upslope::SurfaceWater const > *smartarg1 = 0 ;
+  PyObject *swig_obj[1] ;
+  double result;
+
+  if (!args) SWIG_fail;
+  swig_obj[0] = args;
+  {
+    int newmem = 0;
+    res1 = SWIG_ConvertPtrAndOwn(swig_obj[0], &argp1, SWIGTYPE_p_std__shared_ptrT_cmf__upslope__SurfaceWater_t, 0 |  0 , &newmem);
+    if (!SWIG_IsOK(res1)) {
+      SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "SurfaceWater_get_coverage" "', argument " "1"" of type '" "cmf::upslope::SurfaceWater const *""'"); 
+    }
+    if (newmem & SWIG_CAST_NEW_MEMORY) {
+      tempshared1 = *reinterpret_cast< std::shared_ptr< const cmf::upslope::SurfaceWater > * >(argp1);
+      delete reinterpret_cast< std::shared_ptr< const cmf::upslope::SurfaceWater > * >(argp1);
+      arg1 = const_cast< cmf::upslope::SurfaceWater * >(tempshared1.get());
+    } else {
+      smartarg1 = reinterpret_cast< std::shared_ptr< const cmf::upslope::SurfaceWater > * >(argp1);
+      arg1 = const_cast< cmf::upslope::SurfaceWater * >((smartarg1 ? smartarg1->get() : 0));
+    }
+  }
+  {
+    try {
+      result = (double)((cmf::upslope::SurfaceWater const *)arg1)->get_coverage();
+    } catch (const std::out_of_range& e) {
+      SWIG_exception(SWIG_IndexError, e.what());    
+    } catch (const std::exception& e) {
+      SWIG_exception(SWIG_RuntimeError, e.what());
+    }
+
+  }
+  resultobj = SWIG_From_double(static_cast< double >(result));
+  return resultobj;
+fail:
+  return NULL;
+}
+
+
 SWIGINTERN PyObject *_wrap_SurfaceWater_get_cell(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
   PyObject *resultobj = 0;
   cmf::upslope::SurfaceWater *arg1 = (cmf::upslope::SurfaceWater *) 0 ;
@@ -79196,7 +79189,7 @@ static PyMethodDef SwigMethods[] = {
 		"a cmf::upslope::SurfaceWater \n"
 		""},
 	 { (char *)"Cell_surfacewater_as_storage", (PyCFunction)_wrap_Cell_surfacewater_as_storage, METH_O, (char *)"\n"
-		"Cell_surfacewater_as_storage(Cell self)\n"
+		"Cell_surfacewater_as_storage(Cell self) -> cmf::upslope::surfacewater_ptr\n"
 		"\n"
 		"void surfacewater_as_storage()\n"
 		"\n"
@@ -79251,8 +79244,6 @@ static PyMethodDef SwigMethods[] = {
 		"real albedo()\n"
 		"const \n"
 		""},
-	 { (char *)"Cell_surface_amplitude_set", _wrap_Cell_surface_amplitude_set, METH_VARARGS, (char *)"Cell_surface_amplitude_set(Cell self, real surface_amplitude)"},
-	 { (char *)"Cell_surface_amplitude_get", (PyCFunction)_wrap_Cell_surface_amplitude_get, METH_O, (char *)"Cell_surface_amplitude_get(Cell self) -> real"},
 	 { (char *)"Cell_surface_water_coverage", (PyCFunction)_wrap_Cell_surface_water_coverage, METH_O, (char *)"\n"
 		"Cell_surface_water_coverage(Cell self) -> real\n"
 		"\n"
@@ -80871,6 +80862,7 @@ static PyMethodDef SwigMethods[] = {
 		"Gets the height function (a cmf::river::Prism) for further reference.\n"
 		"\n"
 		""},
+	 { (char *)"SurfaceWater_get_coverage", (PyCFunction)_wrap_SurfaceWater_get_coverage, METH_O, (char *)"SurfaceWater_get_coverage(SurfaceWater self) -> double"},
 	 { (char *)"SurfaceWater_get_cell", (PyCFunction)_wrap_SurfaceWater_get_cell, METH_O, (char *)"\n"
 		"SurfaceWater_get_cell(SurfaceWater self) -> Cell\n"
 		"\n"
@@ -83001,7 +82993,7 @@ static swig_type_info _swigt__p_std__shared_ptrT_cmf__river__Reach_t = {"_p_std_
 static swig_type_info _swigt__p_std__shared_ptrT_cmf__upslope__ET__ShuttleworthWallace_t = {"_p_std__shared_ptrT_cmf__upslope__ET__ShuttleworthWallace_t", "std::shared_ptr< cmf::upslope::ET::ShuttleworthWallace > *", 0, 0, (void*)0, 0};
 static swig_type_info _swigt__p_std__shared_ptrT_cmf__upslope__MacroPore_t = {"_p_std__shared_ptrT_cmf__upslope__MacroPore_t", "std::shared_ptr< cmf::upslope::MacroPore > *", 0, 0, (void*)0, 0};
 static swig_type_info _swigt__p_std__shared_ptrT_cmf__upslope__SoilLayer_t = {"_p_std__shared_ptrT_cmf__upslope__SoilLayer_t", "cmf::upslope::SoilLayer::ptr *|std::shared_ptr< cmf::upslope::SoilLayer > *", 0, 0, (void*)0, 0};
-static swig_type_info _swigt__p_std__shared_ptrT_cmf__upslope__SurfaceWater_t = {"_p_std__shared_ptrT_cmf__upslope__SurfaceWater_t", "std::shared_ptr< cmf::upslope::SurfaceWater > *", 0, 0, (void*)0, 0};
+static swig_type_info _swigt__p_std__shared_ptrT_cmf__upslope__SurfaceWater_t = {"_p_std__shared_ptrT_cmf__upslope__SurfaceWater_t", "std::shared_ptr< cmf::upslope::SurfaceWater > *|cmf::upslope::surfacewater_ptr *", 0, 0, (void*)0, 0};
 static swig_type_info _swigt__p_std__shared_ptrT_cmf__upslope__aquifer_t = {"_p_std__shared_ptrT_cmf__upslope__aquifer_t", "std::shared_ptr< cmf::upslope::aquifer > *", 0, 0, (void*)0, 0};
 static swig_type_info _swigt__p_std__shared_ptrT_cmf__water__DirichletBoundary_t = {"_p_std__shared_ptrT_cmf__water__DirichletBoundary_t", "std::shared_ptr< cmf::water::DirichletBoundary > *", 0, 0, (void*)0, 0};
 static swig_type_info _swigt__p_std__shared_ptrT_cmf__water__NeumannBoundary_t = {"_p_std__shared_ptrT_cmf__water__NeumannBoundary_t", "cmf::water::NeumannBoundary::ptr *|std::shared_ptr< cmf::water::NeumannBoundary > *", 0, 0, (void*)0, 0};

cmf/cmf_core_src/upslope/cell.cpp

@@ -44,7 +44,7 @@ Cell::~Cell()
 Cell::Cell( double _x,double _y,double _z,double area,cmf::project& _project/*=0*/ ) 
 	: x(_x),y(_y),z(_z),m_Area(area),m_project(_project),
 	  m_SurfaceWater(new cmf::water::DirichletBoundary(_project,_z)),	Id(0),
-	  m_meteo(new cmf::atmosphere::ConstantMeteorology), Tground(-300), surface_amplitude(0.01)
+	  m_meteo(new cmf::atmosphere::ConstantMeteorology), Tground(-300)
 {
 	Id = _project.get_cells().size();
 	std::stringstream sstr;
@@ -247,7 +247,7 @@ cmf::water::flux_node::ptr Cell::get_transpiration()
 	return m_Transpiration;
 }
 
-void Cell::surfacewater_as_storage()
+surfacewater_ptr Cell::surfacewater_as_storage()
 {
 	using namespace cmf::river;
 	if (m_SurfaceWaterStorage==0) 
@@ -257,6 +257,7 @@ void Cell::surfacewater_as_storage()
 		m_storages.push_back(m_SurfaceWaterStorage);
 		m_SurfaceWater.reset();
 	}
+	return m_SurfaceWaterStorage;
 }
 
 std::string Cell::to_string() const
@@ -330,6 +331,22 @@ real Cell::heat_flux( cmf::math::Time t) const
 	return Rn + Qs + Ql;
 }
 
+/// Return the fraction of wet leaves in the canopy if a canopy water storage exists.
+/// If no canopy storage is present, it returns 0.0 (=empty). 
+/// The fraction of wet leaves are calculated as the linear filling of the canopy storage.
+
+real cmf::upslope::Cell::leave_wetness() const {
+
+	if (m_Canopy) {
+		// calculate canopy capacity in m3
+		real canopy_capacity = vegetation.LAI * vegetation.CanopyCapacityPerLAI * m_Area * 1e-3;
+		return minmax(m_Canopy->get_volume() / canopy_capacity, 0.0, 1.0);
+	}
+	else {
+		return 0.0;
+	}
+}
+
 real Cell::albedo() const
 {
 	double 
@@ -351,7 +368,7 @@ real Cell::snow_coverage() const
 real Cell::surface_water_coverage() const
 {
 	if (m_SurfaceWaterStorage)
-		return piecewise_linear(m_SurfaceWaterStorage->get_volume()/get_area(),0,surface_amplitude);
+		return m_SurfaceWaterStorage->get_coverage();
 	else
 		return 0.0;
 }

cmf/cmf_core_src/upslope/cell.h

@@ -49,6 +49,9 @@ namespace cmf {
 		typedef void (*connectorfunction)(cmf::upslope::Cell&,cmf::upslope::Cell&,ptrdiff_t);
 		typedef void (*internal_connector)(cmf::upslope::Cell&);
 		typedef std::shared_ptr<SoilLayer> layer_ptr;
+		class SurfaceWater;
+		typedef std::shared_ptr<SurfaceWater> surfacewater_ptr;
+
 		/// A helper class to connect cells with flux_connection objects. This is generated by flux_connection classes, intended to connect cells 
 		class CellConnector
 		{
@@ -117,8 +120,8 @@ namespace cmf {
 			nodemap m_cellboundaries;
 			cmf::water::WaterStorage::ptr
 				m_Canopy,
-				m_Snow,
-				m_SurfaceWaterStorage;
+				m_Snow;
+			cmf::upslope::surfacewater_ptr m_SurfaceWaterStorage;
 			cmf::water::flux_node::ptr m_SurfaceWater;
 			cmf::atmosphere::RainSource::ptr m_rainfall;
 			cmf::water::flux_node::ptr m_Evaporation, m_Transpiration;
@@ -207,7 +210,7 @@ namespace cmf {
 			cmf::water::flux_node::ptr get_surfacewater();
 
 			/// Makes the surfacewater of this cell a cmf::upslope::SurfaceWater storage.
-			void surfacewater_as_storage();
+			surfacewater_ptr surfacewater_as_storage();
 
 			/// Adds a new storage to the cell
 			///
@@ -235,7 +238,6 @@ namespace cmf {
 			cmf::water::WaterStorage::ptr get_snow() const;
 			real snow_coverage() const;
 			real albedo() const;
-			real surface_amplitude;
 			/// Returns the coverage of the surface water.
 			///
 			/// The covered fraction (0..1) is simply modelled as a piecewise linear
@@ -253,17 +255,8 @@ namespace cmf {
 			/// Return the fraction of wet leaves in the canopy if a canopy water storage exists.
 			/// If no canopy storage is present, it returns 0.0 (=empty). 
 			/// The fraction of wet leaves are calculated as the linear filling of the canopy storage.
-			real leave_wetness() const {
-
-				if (m_Canopy) {
-					// calculate canopy capacity in m3
-					real canopy_capacity = vegetation.LAI * vegetation.CanopyCapacityPerLAI * m_Area * 1e-3;
-					return minmax(m_Canopy->get_volume() / canopy_capacity, 0.0, 1.0);
-				}
-				else {
-					return 0.0;
-				}
-			}
+			real leave_wetness() const;
+
 			ptrdiff_t Id;
 			cmf::bytestring get_WKB() const {
 				return m_WKB;

cmf/cmf_core_src/upslope/surfacewater.h

@@ -54,6 +54,11 @@ namespace cmf {
 				throw std::runtime_error("Cannot set the height function of a surface water. Change the parameters");
 			}
 
+			/// Get surface coverage as a function of the actual volume
+			double get_coverage() const {
+				return this->m_height_function->A(this->get_volume()) / this->get_cell().get_area();
+			}
+
 			/// get Manning roughness (n) of the surface
 			///
 			/// From Python use this as a property:

cmf/cmf_core_src/upslope/vegetation/ShuttleworthWallace.cpp

@@ -876,11 +876,8 @@ void cmf::upslope::ET::ShuttleworthWallace::refresh( cmf::math::Time t )
 	// Calculate potential transpiration for wet leave case
 	SWPE(RN,RNground,w.e_s-w.e_a,RAA,RAC,RAS,0,RSS,DELTA,PIR,GER);
 
-	// Get wetness of canopy
-	double VOL_IN_CANOPY = cell.get_canopy() ? cell.m3_to_mm(cell.get_canopy()->get_volume()): 0.0;
-	double MAX_VOL_IN_CANOPY = LAI * cell.vegetation.CanopyCapacityPerLAI;
 	// Calculate actual interception rate from wet leave case and canopy wetness
-	double wet_canopy_fraction = piecewise_linear(VOL_IN_CANOPY, 0, MAX_VOL_IN_CANOPY);
+	double wet_canopy_fraction = cell.leave_wetness();
 	AIR =  wet_canopy_fraction * PIR;
 	// Already for leave evaporation used energy cannot be used for transpiration
 	PTR *= 1 - wet_canopy_fraction;