Index: /issm/trunk-jpl/src/m/sampling.py
===================================================================
--- /issm/trunk-jpl/src/m/sampling.py	(revision 26012)
+++ /issm/trunk-jpl/src/m/sampling.py	(revision 26012)
@@ -0,0 +1,106 @@
+import numpy as np
+
+from checkfield import checkfield
+from fielddisplay import fielddisplay
+from project3d import project3d
+from WriteData import WriteData
+
+
+class sampling(object):
+    """SAMPLING class definition
+
+    Usage:
+        sampling = sampling()
+    """
+
+    def __init__(self):  # {{{
+        self.kappa = float('NaN')
+        self.tau = 0
+        self.beta = float('NaN')
+        self.hydrostatic_adjustment = 0
+        self.phi = 0
+        self.alpha = 0
+        self.robin = 0
+        self.seed = 0
+        self.requested_outputs = []
+
+        # Set defaults
+        self.setdefaultparameters()
+
+    #}}}
+
+    def __repr__(self):  # {{{
+        s = '   Sampling parameters::\n'
+        s += '{}\n'.format(fielddisplay(self,'kappa','coefficient of the identity operator'))
+        s += '{}\n'.format(fielddisplay(self,'tau','scaling coefficient of the solution (default 1.0)'))
+        s += '{}\n'.format(fielddisplay(self,'alpha','exponent in PDE operator, (default 2.0, BiLaplacian covariance operator)'))
+        s += '{}\n'.format(disp(sprintf('\n      %s','Parameters of Robin boundary conditions nabla () \cdot normvec + beta ():')))
+        s += '{}\n'.format(fielddisplay(self,'beta','Coefficient in Robin boundary conditions (to be defined for robin = 1)'))
+        s += '{}\n'.format(fielddisplay(self,'phi','Temporal correlation factor (|phi|<1 for stationary process, phi = 1 for random walk process) (default 0)'))
+        s += '{}\n'.format(fielddisplay(self,'seed','Seed for pseudorandom number generator (given seed if >=0 and random seed if <0) (default -1)'))
+        s += '{}\n'.format(fielddisplay(self, 'requested_outputs', 'additional outputs requested (not implemented yet)'))
+        return s
+    #}}}
+
+    def defaultoutputs(self, md):  # {{{
+        return []
+
+    #}}}
+    def setdefaultparameters(self):  # {{{
+        # Scaling coefficient
+        self.tau = 1
+        # Apply Robin boundary conditions
+        self.robin = 0
+        # Temporal correlation factor
+        self.phi = 0
+        # Exponent in fraction SPDE (default=2, biLaplacian covariance operator)
+        self.alpha=2; 
+        # Seed for pseudorandom number generator (default -1 for random seed)
+        self.alpha=-1; 
+        # Default output
+        self.requested_outputs = ['default']
+        return self
+    #}}}
+
+    def checkconsistency(self, md, solution, analyses):  # {{{
+        # Early return
+        if (SamplingAnalysis' not in analyses):
+            return md
+
+        md = checkfield(md,'fieldname','sampling.kappa','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1],'>',0);
+        md = checkfield(md,'fieldname','sampling.tau','NaN',1,'Inf',1,'numel',1,'>',0);
+        md = checkfield(md,'fieldname','sampling.beta','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1],'>',0);
+        md = checkfield(md,'fieldname','sampling.phi','NaN',1,'Inf',1,'numel',1,'>=',0);
+        md = checkfield(md,'fieldname','sampling.alpha','NaN',1,'Inf',1,'numel',1,'>',0);
+        md = checkfield(md,'fieldname','sampling.robin','numel',1,'values',[0 1]);
+        md = checkfield(md,'fieldname','sampling.seed','NaN',1,'Inf',1,'numel',1);
+        md = checkfield(md,'fieldname','sampling.requested_outputs','stringrow',1);
+
+        return md
+    # }}}
+
+    def marshall(self, prefix, md, fid):  # {{{
+        WriteData(fid,prefix,'object',self,'fieldname','kappa','format','DoubleMat','mattype',1);
+        WriteData(fid,prefix,'object',self,'fieldname','tau','format','Double');
+        WriteData(fid,prefix,'object',self,'fieldname','beta','format','DoubleMat','mattype',1);
+        WriteData(fid,prefix,'object',self,'fieldname','phi','format','Double');
+        WriteData(fid,prefix,'object',self,'fieldname','alpha','format','Integer');
+        WriteData(fid,prefix,'object',self,'fieldname','robin','format','Boolean');
+        WriteData(fid,prefix,'object',self,'fieldname','seed','format','Integer');
+
+        # Process requested outputs
+        outputs = self.requested_outputs
+        indices = [i for i, x in enumerate(outputs) if x == 'default']
+        if len(indices) > 0:
+            outputscopy = outputs[0:max(0, indices[0] - 1)] + self.defaultoutputs(md) + outputs[indices[0] + 1:]
+            outputs = outputscopy
+        WriteData(fid, prefix, 'data', outputs, 'name', 'md.masstransport.requested_outputs', 'format', 'StringArray')
+
+        # Process requested outputs
+		outputs = self.requested_outputs;
+		indices = [i for i, x in enumerate(outputs) if x == 'default']
+		   if len(indices) > 0:
+            outputscopy = outputs[0:max(0, indices[0] - 1)] + self.defaultoutputs(md) + outputs[indices[0] + 1:]
+            outputs = outputscopy
+		WriteData(fid,prefix,'data',outputs,'name','md.sampling.requested_outputs','format','StringArray');
+    # }}}