#35: Removed GIR. RSS now also function of surfacewater

cmf/__init__.py

@@ -22,7 +22,7 @@
 from .stopwatch import StopWatch
 
 
-__version__ = '1.3.1a0.i33_doc_struct'
+__version__ = '1.3.1a0.i35_repair_GIR'
 
 from .cmf_core import connect_cells_with_flux as __ccwf
 

cmf/cmf_core.py

@@ -13299,7 +13299,6 @@ def refresh(self, *args):
     GER = _swig_property(_cmf_core.ShuttleworthWallace_GER_get, _cmf_core.ShuttleworthWallace_GER_set)
     PIR = _swig_property(_cmf_core.ShuttleworthWallace_PIR_get, _cmf_core.ShuttleworthWallace_PIR_set)
     AIR = _swig_property(_cmf_core.ShuttleworthWallace_AIR_get, _cmf_core.ShuttleworthWallace_AIR_set)
-    GIR = _swig_property(_cmf_core.ShuttleworthWallace_GIR_get, _cmf_core.ShuttleworthWallace_GIR_set)
     ATR_sum = _swig_property(_cmf_core.ShuttleworthWallace_ATR_sum_get, _cmf_core.ShuttleworthWallace_ATR_sum_set)
     ATR = _swig_property(_cmf_core.ShuttleworthWallace_ATR_get, _cmf_core.ShuttleworthWallace_ATR_set)
     KSNVP = _swig_property(_cmf_core.ShuttleworthWallace_KSNVP_get, _cmf_core.ShuttleworthWallace_KSNVP_set)

cmf/cmf_core_src/cmf_wrap.cpp

@@ -73408,58 +73408,6 @@ SWIGINTERN PyObject *_wrap_ShuttleworthWallace_AIR_get(PyObject *SWIGUNUSEDPARM(
 }
 
 
-SWIGINTERN PyObject *_wrap_ShuttleworthWallace_GIR_set(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
-  PyObject *resultobj = 0;
-  cmf::upslope::ET::ShuttleworthWallace *arg1 = (cmf::upslope::ET::ShuttleworthWallace *) 0 ;
-  double arg2 ;
-  void *argp1 = 0 ;
-  int res1 = 0 ;
-  double val2 ;
-  int ecode2 = 0 ;
-  PyObject *swig_obj[2] ;
-
-  if (!SWIG_Python_UnpackTuple(args,"ShuttleworthWallace_GIR_set",2,2,swig_obj)) SWIG_fail;
-  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_cmf__upslope__ET__ShuttleworthWallace, 0 |  0 );
-  if (!SWIG_IsOK(res1)) {
-    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "ShuttleworthWallace_GIR_set" "', argument " "1"" of type '" "cmf::upslope::ET::ShuttleworthWallace *""'"); 
-  }
-  arg1 = reinterpret_cast< cmf::upslope::ET::ShuttleworthWallace * >(argp1);
-  ecode2 = SWIG_AsVal_double(swig_obj[1], &val2);
-  if (!SWIG_IsOK(ecode2)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "ShuttleworthWallace_GIR_set" "', argument " "2"" of type '" "double""'");
-  } 
-  arg2 = static_cast< double >(val2);
-  if (arg1) (arg1)->GIR = arg2;
-  resultobj = SWIG_Py_Void();
-  return resultobj;
-fail:
-  return NULL;
-}
-
-
-SWIGINTERN PyObject *_wrap_ShuttleworthWallace_GIR_get(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
-  PyObject *resultobj = 0;
-  cmf::upslope::ET::ShuttleworthWallace *arg1 = (cmf::upslope::ET::ShuttleworthWallace *) 0 ;
-  void *argp1 = 0 ;
-  int res1 = 0 ;
-  PyObject *swig_obj[1] ;
-  double result;
-
-  if (!args) SWIG_fail;
-  swig_obj[0] = args;
-  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_cmf__upslope__ET__ShuttleworthWallace, 0 |  0 );
-  if (!SWIG_IsOK(res1)) {
-    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "ShuttleworthWallace_GIR_get" "', argument " "1"" of type '" "cmf::upslope::ET::ShuttleworthWallace *""'"); 
-  }
-  arg1 = reinterpret_cast< cmf::upslope::ET::ShuttleworthWallace * >(argp1);
-  result = (double) ((arg1)->GIR);
-  resultobj = SWIG_From_double(static_cast< double >(result));
-  return resultobj;
-fail:
-  return NULL;
-}
-
-
 SWIGINTERN PyObject *_wrap_ShuttleworthWallace_ATR_sum_set(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
   PyObject *resultobj = 0;
   cmf::upslope::ET::ShuttleworthWallace *arg1 = (cmf::upslope::ET::ShuttleworthWallace *) 0 ;
@@ -81638,8 +81586,6 @@ static PyMethodDef SwigMethods[] = {
 	 { (char *)"ShuttleworthWallace_PIR_get", (PyCFunction)_wrap_ShuttleworthWallace_PIR_get, METH_O, (char *)"ShuttleworthWallace_PIR_get(ShuttleworthWallace self) -> double"},
 	 { (char *)"ShuttleworthWallace_AIR_set", _wrap_ShuttleworthWallace_AIR_set, METH_VARARGS, (char *)"ShuttleworthWallace_AIR_set(ShuttleworthWallace self, double AIR)"},
 	 { (char *)"ShuttleworthWallace_AIR_get", (PyCFunction)_wrap_ShuttleworthWallace_AIR_get, METH_O, (char *)"ShuttleworthWallace_AIR_get(ShuttleworthWallace self) -> double"},
-	 { (char *)"ShuttleworthWallace_GIR_set", _wrap_ShuttleworthWallace_GIR_set, METH_VARARGS, (char *)"ShuttleworthWallace_GIR_set(ShuttleworthWallace self, double GIR)"},
-	 { (char *)"ShuttleworthWallace_GIR_get", (PyCFunction)_wrap_ShuttleworthWallace_GIR_get, METH_O, (char *)"ShuttleworthWallace_GIR_get(ShuttleworthWallace self) -> double"},
 	 { (char *)"ShuttleworthWallace_ATR_sum_set", _wrap_ShuttleworthWallace_ATR_sum_set, METH_VARARGS, (char *)"ShuttleworthWallace_ATR_sum_set(ShuttleworthWallace self, double ATR_sum)"},
 	 { (char *)"ShuttleworthWallace_ATR_sum_get", (PyCFunction)_wrap_ShuttleworthWallace_ATR_sum_get, METH_O, (char *)"ShuttleworthWallace_ATR_sum_get(ShuttleworthWallace self) -> double"},
 	 { (char *)"ShuttleworthWallace_ATR_set", _wrap_ShuttleworthWallace_ATR_set, METH_VARARGS, (char *)"ShuttleworthWallace_ATR_set(ShuttleworthWallace self, cmf::math::num_array ATR)"},

cmf/cmf_core_src/upslope/vegetation/ShuttleworthWallace.cpp

@@ -756,7 +756,7 @@ void SNOVAP (double TSNOW, double TA, double EA, double UA, double ZA, double HE
 	PSNVP = KSNVP * PSNVP;
 }
 cmf::upslope::ET::ShuttleworthWallace::ShuttleworthWallace( cmf::upslope::Cell& _cell, bool _allow_dew) 
-:	PTR(0.0), ATR_sum(0.0), ATR(), GER(0.0),PIR(0.0),GIR(0.0), PSNVP(0.0), ASNVP(0.0), cell(_cell),KSNVP(1.0),AIR(0.0),
+:	PTR(0.0), ATR_sum(0.0), ATR(), GER(0.0),PIR(0.0), PSNVP(0.0), ASNVP(0.0), cell(_cell),KSNVP(1.0),AIR(0.0),
 	RAA(0.0),RAC(0.0),RAS(0.0),RSS(0.0),RSC(0.0), refresh_counter(0), allow_dew(_allow_dew)
 {
 	KSNVP = piecewise_linear(cell.vegetation.Height,1,5,1.0,0.3);
@@ -842,34 +842,51 @@ void cmf::upslope::ET::ShuttleworthWallace::refresh( cmf::math::Time t )
 	double albedo = cell.snow_coverage() * 0.9 + (1-cell.snow_coverage()) * v.albedo;
 	double RN = w.Rn(albedo) / WTOMJ;
 	double RNground = RN * exp(-v.CanopyPARExtinction * (LAI+SAI));
-	RSS = FRSS(RSSa ,RSSb, RSSa_pot, cell.get_layer(0)->get_matrix_potential());
+
+	// Get atmospheric resitances RAA, RAC, RAS
 	SWGRA(UA,ZA,h,Z0,DISP,Z0C,DISPC,Z0GS,v.LeafWidth,RHOTP,NN,LAI,SAI,RAA,RAC,RAS);
-	if (w.Rs>0) {
+
+	// Get canopy surface resistance
+	if (w.Rs>0) { // by day
 		SRSC(w.Rs / WTOMJ, w.T, w.e_s-w.e_a, LAI, SAI, GLMIN, GLMAX,
 			R5,CVPD,RM,v.CanopyPARExtinction,
 			T_f[0],T_f[1],T_f[2],T_f[3],
 			RSC);
-	} else {
+	} else { //by night
 		RSC = 1 / (GLMIN * LAI);
 	}
+
+	// Get soil surface resitance neglecting surface water
+	RSS = FRSS(RSSa, RSSb, RSSa_pot, cell.get_layer(0)->get_matrix_potential());
+	// Lower RSS for surfacewater (zero 
+	RSS *= 1 - cell.surface_water_coverage();
+
+	// Get saturated vapor pressure factor
 	double DELTA,ES;
 	ESAT(w.T,ES,DELTA);
 
+	// Calculate potential transpiration for dry leave case
 	SWPE(RN,RNground,w.e_s-w.e_a,RAA,RAC,RAS,RSC,RSS,DELTA,PTR,GER);
-	// Shuttleworth-Wallace potential interception and ground evap. rates
-	// RSC = 0, RSS not changed
-	SWPE(RN,RNground,w.e_s-w.e_a,RAA,RAC,RAS,0,RSS,DELTA,PIR,GIR);
+	
+	// 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
 	AIR =  piecewise_linear(VOL_IN_CANOPY,0,MAX_VOL_IN_CANOPY,0,PIR);
+	// Already for leave evaporation used energy cannot be used for transpiration
 	PTR -= AIR;
+	// If energy for transpiration is left (dry or partial dry leaves)
 	if (PTR>0.001)	{
+		// Calculate actual transpiration
 		TBYLAYER(1,PTR,DISPC,ALPHA,KK,RROOTI,RXYLEM,PSITI,lc, RHOWG * pF_to_waterhead(4.2),1,ATR_sum,ATR);
 		if (ATR_sum < PTR)
+			// Update ground evaporation
 			SWGE(RN,RNground,w.e_s-w.e_a,RAA,RAS,RSS,DELTA,ATR_sum + AIR,GER);
 
-	} else {
+	} else { // wet leaves no energy for transpiration
 		PTR = 0.0;
 		ATR = 0.0;
 		ATR_sum=0.0;
@@ -880,7 +897,6 @@ void cmf::upslope::ET::ShuttleworthWallace::refresh( cmf::math::Time t )
 		PTR = maximum(PTR,0.0);
 		GER = maximum(GER,0.0);
 		PIR = maximum(PIR,0.0);
-		GIR = maximum(GIR,0.0);
 	}
 
 }
@@ -921,7 +937,7 @@ double cmf::upslope::ET::ShuttleworthWallace::evap_from_openwater( cmf::river::O
 		res += (1-cmf::upslope::connections::snowfraction(w.T)) * GER * 1e-3 * cell.get_area();
 	} 
 	// Get evaporation from open water, if open water is not empty
-	res += GIR * ows->get_height_function().A(ows->get_volume()) * 1e-3 * (1-ows->is_empty());
+	res += GER * cell.surface_water_coverage();
 	return res;
 }
 

cmf/cmf_core_src/upslope/vegetation/ShuttleworthWallace.h

@@ -93,15 +93,13 @@ namespace cmf {
 				double PSNVP;
 				/// Actual snow evaporation rate in mm/day
 				double ASNVP;
-				/// Drought stressed soil evaporation rate in mm/day (actual rate may be lower due to snow or surface water cover)
+				/// Ground evaporation rate in mm/day (either from soil or from surfacewater)
 				double GER; 
 				/// Potential leaf evaporation rate in mm/day
 				double PIR;
 				/// Actual leaf evaporation rate in mm/day
 				double AIR;
 				/// potential surface water evaporation rate in mm/day, actual rate depends on surfacewater coverage
-				double GIR; 
-				/// actual transpiration rate in mm/day
 				double ATR_sum;
 				/// actual transpiration rate per layer in mm/day
 				cmf::math::num_array ATR;

demo/surface-sw-et.py

@@ -4,11 +4,13 @@
 """
 import cmf
 
+
 class Model:
+
     def __init__(self):
         self.p = cmf.project()
         self.c = self.p.NewCell(0, 0, 0, 1000, with_surfacewater=True)
-        l = self.c.add_layer(1.0)
+        self.c.add_layer(1.0)
         self.c: cmf.Cell
         cmf.RutterInterception.use_for_cell(self.c)
         self.et = cmf.ShuttleworthWallace.use_for_cell(self.c)
@@ -17,7 +19,7 @@ def __init__(self):
     def __call__(self, days=7):
         self.c.vegetation.LAI = 7
         self.c.vegetation.height = 10
-        self.c.vegetation.CanopyClosure = 0.5
+        self.c.vegetation.CanopyClosure = 1
         l = self.c.layers[0]
         self.c.surfacewater.depth = 0.05
         l.soil.Ksat = 0.02
@@ -30,8 +32,9 @@ def __call__(self, days=7):
         solver = cmf.CVodeIntegrator(self.p, 1e-9)
         vol = []
         flux = []
+        resistance = []
+        mpot = []
         self.et.refresh(cmf.Time())
-        print(f'GER={self.et.GER:0.3} mm, GIR={self.et.GIR:0.3} mm')
         for t in solver.run(cmf.Time(), cmf.day * days, cmf.min*6):
             #print(f'{t}: {self.c.surfacewater.depth:0.3f}m')
             vol.append((
@@ -46,19 +49,28 @@ def __call__(self, days=7):
                 self.c.layers[0].flux_to(self.c.transpiration, t),
                 self.c.surfacewater.flux_to(self.c.layers[0], t),
             ))
-        return vol, flux
+            resistance.append((
+                self.et.RAA, self.et.RAC, self.et.RAS,
+                self.et.RSC, self.et.RSS
+            ))
+            mpot.append(self.c.surface_water_coverage())
+        return vol, flux, resistance, mpot
 
 if __name__ == '__main__':
+    print(cmf.__version__)
     m = Model()
-    vol, flux = m(7)
+    vol, flux, resistance, mpot = m(10)
     from matplotlib import pylab as plt
-    fig, ax = plt.subplots(2, 1, True)
+    fig, ax = plt.subplots(3, 1, sharex='all')
     plt.sca(ax[0])
     plt.plot(vol)
-    plt.legend(['Canopy', 'Surfacewater', 'Layer'])
+    plt.legend(['Canopy', 'Surfacewater', 'Layer'], loc=0)
     plt.sca(ax[1])
     plt.plot(flux)
-    plt.legend(['E_Canopy', 'E_Surfacewater', 'E_Layer', 'T_Layer', 'Infiltration'])
+    plt.legend(['E_Canopy', 'E_Surfacewater', 'E_Layer', 'T_Layer', 'Infiltration'], loc=0)
+    plt.sca(ax[2])
+    plt.plot(resistance)
+    plt.legend('RAA RAC RAS RSC RSS'.split(), loc=0)
     plt.show()