Index: /issm/trunk-jpl/src/m/classes/model.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/model.m	(revision 25995)
+++ /issm/trunk-jpl/src/m/classes/model.m	(revision 25996)
@@ -39,10 +39,11 @@
 		steadystate      = 0;
 		transient        = 0;
-		levelset			  = 0;
+		levelset	     = 0;
 		calving          = 0;
 		frontalforcings  = 0;
-		love			     = 0;
-		gia				  = 0;
+		love			 = 0;
+		gia				 = 0;
 		esa              = 0;
+        sampling         = 0;
 
 		autodiff         = 0;
@@ -187,5 +188,7 @@
 					end
 				end
-			end
+            end
+            %2021 February 17
+			if isa(md.sampling,'double'); md.sampling=sampling(); end
 		end% }}}
 	end
@@ -241,4 +244,5 @@
 			md.love             = fourierlove();
 			md.esa              = esa();
+            md.sampling         = sampling();
 			md.autodiff         = autodiff();
 			md.inversion        = inversion();
@@ -1505,4 +1509,5 @@
 			disp(sprintf('%19s: %-22s -- %s','esa'             ,['[1x1 ' class(self.esa) ']'],'parameters for elastic adjustment solution'));
 			disp(sprintf('%19s: %-22s -- %s','love'            ,['[1x1 ' class(self.love) ']'],'parameters for love solution'));
+            disp(sprintf('%19s: %-22s -- %s','sampling'        ,['[1x1 ' class(self.sampling) ']'],'parameters for stochastic sampler'));
 			disp(sprintf('%19s: %-22s -- %s','autodiff'        ,['[1x1 ' class(self.autodiff) ']'],'automatic differentiation parameters'));
 			disp(sprintf('%19s: %-22s -- %s','inversion'       ,['[1x1 ' class(self.inversion) ']'],'parameters for inverse methods'));
Index: /issm/trunk-jpl/src/m/classes/sampling.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/sampling.m	(revision 25996)
+++ /issm/trunk-jpl/src/m/classes/sampling.m	(revision 25996)
@@ -0,0 +1,129 @@
+%SAMPLING class definition
+%
+%   Usage:
+%      sampling=sampling();
+
+classdef sampling
+	properties (SetAccess=public)
+         kappa          = NaN;
+		 tau            = 0;
+		 beta           = NaN;
+         phi            = 0;
+         alpha          = 0;
+         robin          = 0;
+         seed           = 0;
+		 requested_outputs      = {};
+    end
+	methods
+		function self = sampling(varargin) % {{{
+            switch nargin
+				case 0
+					self=setdefaultparameters(self);
+				otherwise
+					error('constructor not supported');
+            end
+		end % }}}
+		function list = defaultoutputs(self,md) % {{{
+
+            list = {};
+
+		end % }}}
+		function self = 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; % Default 
+            
+            %Seed for pseudorandom number generator (default -1 for random seed)
+            self.seed=-1;
+            
+			%default output
+			self.requested_outputs={'default'};
+
+		end % }}}
+        function md = setparameters(self,md,lc,sigma) % {{{
+            
+            nu = self.alpha-1; 
+            KAPPA = sqrt(8*nu)/lc;
+            TAU = sqrt(gamma(nu)/(gamma(self.alpha)*(4*pi)*KAPPA^(2*nu)*sigma^2));
+            md.sampling.kappa = KAPPA*ones(md.mesh.numberofvertices,1);
+            md.sampling.tau = TAU;
+       
+		end % }}}
+		function md = checkconsistency(self,md,solution,analyses) % {{{
+            
+            if ~ismember('SamplingAnalysis',analyses), return; end
+			
+            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);
+
+		end % }}}
+		function disp(self) % {{{
+
+            disp(sprintf('   Sampling parameters:'));
+
+			disp(sprintf('\n      %s','Parameters of PDE operator (kappa^2 I-Laplacian)^(alpha/2)(tau):'));
+            fielddisplay(self,'kappa','coefficient of the identity operator');
+			fielddisplay(self,'tau','scaling coefficient of the solution (default 1.0)');
+			fielddisplay(self,'alpha','exponent in PDE operator, (default 2.0, BiLaplacian covariance operator)');
+          
+			disp(sprintf('\n      %s','Parameters of Robin boundary conditions nabla () \cdot normvec + beta ():'));
+			fielddisplay(self,'robin','Apply Robin boundary conditions (1 if applied and 0 for homogenous Neumann boundary conditions) (default 0)');
+			fielddisplay(self,'beta','Coefficient in Robin boundary conditions (to be defined for robin = 1)');          
+            
+            disp(sprintf('\n      %s','Parameters for first-order autoregressive process (X_t = phi X_{t-1} + noise) (if transient):'));
+			fielddisplay(self,'phi','Temporal correlation factor (|phi|<1 for stationary process, phi = 1 for random walk process) (default 0)');
+            
+			disp(sprintf('\n      %s','Other parameters of stochastic sampler:'));
+            fielddisplay(self,'seed','Seed for pseudorandom number generator (given seed if >=0 and random seed if <0) (default -1)');
+			fielddisplay(self,'requested_outputs','additional outputs requested (not implemented yet)');
+ 
+		end % }}}
+		function 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;
+			pos  = find(ismember(outputs,'default'));
+			if ~isempty(pos),
+				outputs(pos) = [];                         %remove 'default' from outputs
+				outputs      = [outputs defaultoutputs(self,md)]; %add defaults
+			end
+			WriteData(fid,prefix,'data',outputs,'name','md.sampling.requested_outputs','format','StringArray');
+		end % }}}
+		function savemodeljs(self,fid,modelname) % {{{
+            
+            writejsdouble(fid,[modelname '.sampling.kappa'],self.kappa);
+            writejsdouble(fid,[modelname '.sampling.tau'],self.tau);
+            writejsdouble(fid,[modelname '.sampling.beta'],self.beta);
+            writejsdouble(fid,[modelname '.sampling.phi'],self.beta);
+            writejsdouble(fid,[modelname '.sampling.alpha'],self.alpha);
+            writejsdouble(fid,[modelname '.sampling.robin'],self.robin);
+            writejsdouble(fid,[modelname '.sampling.seed'],self.seed);
+			writejscellstring(fid,[modelname '.sampling.requested_outputs'],self.requested_outputs);
+
+		end % }}}
+	end
+end
Index: /issm/trunk-jpl/src/m/consistency/ismodelselfconsistent.m
===================================================================
--- /issm/trunk-jpl/src/m/consistency/ismodelselfconsistent.m	(revision 25995)
+++ /issm/trunk-jpl/src/m/consistency/ismodelselfconsistent.m	(revision 25996)
@@ -78,4 +78,6 @@
 	elseif strcmp(solutiontype,'DamageEvolutionSolution')
 		analyses={'DamageEvolutionAnalysis'};
+    elseif strcmp(solutiontype,'SamplingSolution')
+		analyses={'SamplingAnalysis'};    
 	else
 		error(' solution type: %s' , solutiontype, ' not supported yet!');