diff --git a/pub/src/curve.cpp b/pub/src/curve.cpp
index e2f158aaf2df9075190e3c4f0e2da632e96f9c22..ad69b2cb23e7edcd2f7e7f1e6afd6aedd38d22fd 100644
--- a/pub/src/curve.cpp
+++ b/pub/src/curve.cpp
@@ -45,8 +45,8 @@ using boost::format;
 namespace NCurveFile
 {
     // Local functions (wouldn't suffice "static" ??):
-    double NumericIntegral(
-        POlc fc, int k, int j, double low, double hig );
+    double NumericIntegral( POlc fc, int k, int j, double low, double hig );
+    double GridSum( POlc fc, int k, int j );
     // Parameters for numeric integration:
     double numint_bound_low=-INFINITY, numint_bound_hig=INFINITY;
     // Tuning parameters:
@@ -1024,6 +1024,7 @@ void NCurveFile::IntegralProperty()
             "   0 | 1 | ...         extract curve parameters\n"
             "   q0 | q1 | q2        extract fit metrics\n"
             "   i                   numeric integral\n"
+            "   s                   integral as sum over data file grid\n"
             "   <list>              comma-separated list thereof\n"
             "   *                   extract all curve parameters\n"
             ;
@@ -1129,6 +1130,8 @@ void NCurveFile::IntegralProperty()
                 else if ( mode=="i" )
                     val = NumericIntegral( fin, fiter.k(), j,
                                            numint_bound_low, numint_bound_hig );
+                else if ( mode=="s" )
+                    val = GridSum( fin, fiter.k(), j );
                 else
                     throw "invalid mode " + mode;
 
@@ -1172,6 +1175,31 @@ double NCurveFile::NumericIntegral(
                             numint_mode, numint_epsabs, numint_epsrel );
 }
 
+//! Numeric integration of one curve, as a sum using the grid from data file kd.
+
+double NCurveFile::GridSum( POlc fc, int k, int j )
+{
+    if( fc->kd==-1 )
+        throw "data reference not set";
+    POld fd = NOlm::mem_get_D(fc->kd);
+    if( !fd )
+        throw "data reference not found";
+    if( j>=fd->nJ() )
+        throw "number of spectra does not match";
+    vector<double> *xc = &(fd->VS(j)->x);
+    int n = xc->size();
+    if( n<2 )
+        throw "not enough points in grid";
+    vector<double> yc;
+    fc->curve_val_vec( &yc, *xc, k, j );
+    double ret = 0;
+    ret  = ((*xc)[1]-(*xc)[0])/2*yc[0];
+    ret += ((*xc)[n-1]-(*xc)[n-2])/2*yc[n-1];
+    for ( uint i=1; i<n-1; ++i )
+        ret += ((*xc)[i+1]-(*xc)[i-1])/2*yc[i];
+    return ret;
+}
+
 
 //***************************************************************************//
 //* Calculate cumulative distribution (needed for SPHERES PST simulation)   *//