Index: issm/trunk-jpl/src/m/classes/slr.m
===================================================================
--- issm/trunk-jpl/src/m/classes/slr.m	(revision 22808)
+++ issm/trunk-jpl/src/m/classes/slr.m	(revision 22955)
@@ -8,4 +8,5 @@
 		deltathickness = NaN;
 		sealevel       = NaN; 
+		spcthickness   = NaN; 
 		maxiter        = 0;
 		reltol         = 0;
@@ -25,5 +26,11 @@
 		ocean_area_scaling = 0;
 		steric_rate    = 0; %rate of ocean expansion from steric effects. 
+		geodetic_run_frequency  = 1; %how many time steps we skip before we run the geodetic part of the solver during transient
+		geodetic       = 0; %compute geodetic SLR? (in addition to steric?)
 		degacc         = 0;
+		loop_increment = 0;
+		horiz = 0;
+		Ngia  = 0;
+		Ugia  = 0;
 		requested_outputs      = {};
 		transitions    = {};
@@ -46,10 +53,12 @@
 		%maximum of non-linear iterations.
 		self.maxiter=5;
+		self.loop_increment=200;
 
 		%computational flags: 
+		self.geodetic=0;
 		self.rigid=1;
 		self.elastic=1;
-		self.rotation=0;
 		self.ocean_area_scaling=0;
+		self.rotation=1;
 
 		%tidal love numbers: 
@@ -73,4 +82,6 @@
 		self.steric_rate=0;
 	
+		%how many time steps we skip before we run SLR solver during transient
+		self.geodetic_run_frequency=1;
 	
 		%output default:
@@ -79,4 +90,7 @@
 		%transitions should be a cell array of vectors: 
 		self.transitions={};
+
+		%horizontal displacement?  (not by default)
+		self.horiz=0;
 		
 		end % }}}
@@ -86,4 +100,5 @@
 			md = checkfield(md,'fieldname','slr.deltathickness','NaN',1,'Inf',1,'size',[md.mesh.numberofelements 1]);
 			md = checkfield(md,'fieldname','slr.sealevel','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
+			md = checkfield(md,'fieldname','slr.spcthickness','Inf',1,'timeseries',1);
 			md = checkfield(md,'fieldname','slr.love_h','NaN',1,'Inf',1);
 			md = checkfield(md,'fieldname','slr.love_k','NaN',1,'Inf',1);
@@ -98,7 +113,12 @@
 			md = checkfield(md,'fieldname','slr.abstol','size',[1 1]);
 			md = checkfield(md,'fieldname','slr.maxiter','size',[1 1],'>=',1);
+			md = checkfield(md,'fieldname','slr.geodetic_run_frequency','size',[1 1],'>=',1);
 			md = checkfield(md,'fieldname','slr.steric_rate','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
 			md = checkfield(md,'fieldname','slr.degacc','size',[1 1],'>=',1e-10);
 			md = checkfield(md,'fieldname','slr.requested_outputs','stringrow',1);
+			md = checkfield(md,'fieldname','slr.loop_increment','NaN',1,'Inf',1,'>=',1);
+			md = checkfield(md,'fieldname','slr.horiz','NaN',1,'Inf',1,'values',[0 1]);
+			md = checkfield(md,'fieldname','slr.Ngia','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
+			md = checkfield(md,'fieldname','slr.Ugia','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
 
 			%check that love numbers are provided at the same level of accuracy: 
@@ -112,6 +132,21 @@
 			[els,vertices]=find(maskpos>0);
 			if length(els),
-				error('slr checkconsistency fail: there are elements with ice loads where some vertices are not on the ice!');
-			end
+				warning('slr checkconsistency fail: there are elements with ice loads where some vertices are not on the ice!');
+			end
+
+			%check that  if geodetic is requested, we are a mesh3dsurface model (planet), or if we are not, 
+			%a coupler to a planet model is provided. 
+			if self.geodetic,
+				if md.transient.iscoupler, 
+					%we are good; 
+				else
+					if strcmpi(class(md.mesh),'mesh3dsurface'),
+						%we are good
+					else
+						error('model is requesting geodetic computations without being a mesh3dsurface, or being coupled to one!');
+					end
+				end
+			end
+
 
 		end % }}}
@@ -124,6 +159,7 @@
 			fielddisplay(self,'deltathickness','thickness change: ice height equivalent [m]');
 			fielddisplay(self,'sealevel','current sea level (prior to computation) [m]');
-			fielddisplay(self,'reltol','sea level rise relative convergence criterion, (NaN: not applied)');
-			fielddisplay(self,'abstol','sea level rise absolute convergence criterion, (default, NaN: not applied)');
+			fielddisplay(self,'spcthickness','thickness constraints (NaN means no constraint) [m]');
+			fielddisplay(self,'reltol','sea level rise relative convergence criterion, (default, NaN: not applied)');
+			fielddisplay(self,'abstol','sea level rise absolute convergence criterion, NaN: not applied');
 			fielddisplay(self,'maxiter','maximum number of nonlinear iterations');
 			fielddisplay(self,'love_h','load Love number for radial displacement');
@@ -136,9 +172,14 @@
 			fielddisplay(self,'polar_moi','polar moment of inertia [kg m^2]');
 			fielddisplay(self,'angular_velocity','mean rotational velocity of earth [per second]'); 
+			fielddisplay(self,'ocean_area_scaling','correction for model representation of ocean area [default: No correction]'); 
+			fielddisplay(self,'steric_rate','rate of steric ocean expansion (in mm/yr)'); 
+			fielddisplay(self,'Ngia','rate of viscous (GIA) geoid expansion (in mm/yr)'); 
+			fielddisplay(self,'Ugia','rate of viscous (GIA) bedrock uplift (in mm/yr)'); 
+			fielddisplay(self,'loop_increment','vector assembly (in the convolution) framentation');
+			fielddisplay(self,'geodetic','compute geodetic SLR? ( in addition to steric?) default 0');
+			fielddisplay(self,'geodetic_run_frequency','how many time steps we skip before we run SLR solver during transient (default: 1)');
 			fielddisplay(self,'rigid','rigid earth graviational potential perturbation');
 			fielddisplay(self,'elastic','elastic earth graviational potential perturbation');
 			fielddisplay(self,'rotation','earth rotational potential perturbation');
-			fielddisplay(self,'ocean_area_scaling','correction for model representation of ocean area (default: No correction)'); 
-			fielddisplay(self,'steric_rate','rate of steric ocean expansion [mm/yr]'); 
 			fielddisplay(self,'degacc','accuracy (default .01 deg) for numerical discretization of the Green''s functions');
 			fielddisplay(self,'transitions','indices into parts of the mesh that will be icecaps');
@@ -149,5 +190,7 @@
 			%WriteData(fid,prefix,'object',self,'fieldname','deltathickness','format','DoubleMat','mattype',2);
 			WriteData(fid,prefix,'object',self,'fieldname','deltathickness','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts);
+			%WriteData(fid,prefix,'object',self,'fieldname','deltathickness','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofelements+1);
 			WriteData(fid,prefix,'object',self,'fieldname','sealevel','mattype',1,'format','DoubleMat','timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
+			WriteData(fid,prefix,'object',self,'fieldname','spcthickness','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
 			WriteData(fid,prefix,'object',self,'fieldname','reltol','format','Double');
 			WriteData(fid,prefix,'object',self,'fieldname','abstol','format','Double');
@@ -166,7 +209,13 @@
 			WriteData(fid,prefix,'object',self,'fieldname','rotation','format','Boolean');
 			WriteData(fid,prefix,'object',self,'fieldname','ocean_area_scaling','format','Boolean');
+			WriteData(fid,prefix,'object',self,'fieldname','geodetic_run_frequency','format','Integer');
 			WriteData(fid,prefix,'object',self,'fieldname','steric_rate','format','DoubleMat','mattype',1,'scale',1e-3/md.constants.yts);
+			WriteData(fid,prefix,'object',self,'fieldname','Ngia','format','DoubleMat','mattype',1,'scale',1e-3/md.constants.yts);
+			WriteData(fid,prefix,'object',self,'fieldname','Ugia','format','DoubleMat','mattype',1,'scale',1e-3/md.constants.yts);
 			WriteData(fid,prefix,'object',self,'fieldname','degacc','format','Double');
 			WriteData(fid,prefix,'object',self,'fieldname','transitions','format','MatArray');
+			WriteData(fid,prefix,'object',self,'fieldname','loop_increment','format','Integer');
+			WriteData(fid,prefix,'object',self,'fieldname','horiz','format','Integer');
+			WriteData(fid,prefix,'object',self,'fieldname','geodetic','format','Integer');
 			
 			%process requested outputs
@@ -184,4 +233,5 @@
 			writejs1Darray(fid,[modelname '.slr.deltathickness'],self.deltathickness);
 			writejs1Darray(fid,[modelname '.slr.sealevel'],self.sealevel);
+			writejs1Darray(fid,[modelname '.slr.spcthickness'],self.spcthickness);
 			writejsdouble(fid,[modelname '.slr.maxiter'],self.maxiter);
 			writejsdouble(fid,[modelname '.slr.reltol'],self.reltol);
@@ -200,4 +250,5 @@
 			writejsdouble(fid,[modelname '.slr.rotation'],self.rotation);
 			writejsdouble(fid,[modelname '.slr.ocean_area_scaling'],self.ocean_area_scaling);
+			writejsdouble(fid,[modelname '.slr.geodetic_run_frequency'],self.geodetic_run_frequency);
 			writejs1Darray(fid,[modelname '.slr.steric_rate'],self.steric_rate);
 			writejsdouble(fid,[modelname '.slr.degacc'],self.degacc);
@@ -205,4 +256,7 @@
 			writejscellarray(fid,[modelname '.slr.transitions'],self.transitions);
 		end % }}}
+		function self = extrude(self,md) % {{{
+			self.sealevel=project3d(md,'vector',self.sealevel,'type','node');
+		end % }}}
 	end
 end
