Index: /issm/trunk-jpl/src/m/classes/slr.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/slr.py	(revision 20158)
+++ /issm/trunk-jpl/src/m/classes/slr.py	(revision 20158)
@@ -0,0 +1,125 @@
+from fielddisplay import fielddisplay
+from MatlabFuncs import *
+from model import *
+from EnumDefinitions import *
+from numpy import *
+from checkfield import checkfield
+from WriteData import WriteData
+
+class slr(object):
+	"""
+	SLR class definition
+	
+		Usage:
+		  slr=slr();
+	"""
+	
+	def __init__(self): # {{{
+		self.deltathickness 	= NaN
+		self.maxiter			= 0
+		self.reltol       		= 0
+		self.bstol         		= 0
+		self.love_h         	= 0 #provided by PREM model()
+		self.love_k         	= 0 #ideam
+		self.rigid         		= 0
+		self.elastic         	= 0
+		self.eustatic         	= 0
+		self.degacc         	= 0
+		self.requested_outputs 	= []
+		self.transitions    	= []
+		
+		#set defaults
+		self.setdefaultparameters()
+		#}}}
+	def __repr__(self): # {{{
+			string='   slr parameters:'
+			string="%s\n%s"%(string,fielddisplay(self,'deltathickness','thickness change (main loading of the slr solution core [m]'))
+			string="%s\n%s"%(string,fielddisplay(self,'reltol','sea level rise relative convergence criterion, (default, NaN: not applied)'))
+			string="%s\n%s"%(string,fielddisplay(self,'abstol','sea level rise absolute convergence criterion, NaN: not applied'))
+			string="%s\n%s"%(string,fielddisplay(self,'maxiter','maximum number of nonlinear iterations'))
+			string="%s\n%s"%(string,fielddisplay(self,'love_h','love load number for radial displacement'))
+			string="%s\n%s"%(string,fielddisplay(self,'love_k','love load number for gravitational potential perturbation'))
+			string="%s\n%s"%(string,fielddisplay(self,'rigid','rigid earth graviational potential perturbation'))
+			string="%s\n%s"%(string,fielddisplay(self,'elastic','elastic earth graviational potential perturbation'))
+			string="%s\n%s"%(string,fielddisplay(self,'eustatic','eustatic sea level rise'))
+			string="%s\n%s"%(string,fielddisplay(self,'degacc','accuracy (default .01 deg) for numerical discretization of the Green''s functions'))
+			string="%s\n%s"%(string,fielddisplay(self,'transitions','indices into parts of the mesh that will be icecaps'))
+			string="%s\n%s"%(string,fielddisplay(self,'requested_outputs','additional outputs requested'))
+
+			return string
+		# }}}
+	def setdefaultparameters(self): # {{{
+		
+		#Convergence criterion: absolute, relative and residual
+		self.reltol=NaN #default
+		self.abstol=0.001 #1 mm of sea level rise
+
+		#maximum of non-linear iterations.
+		self.maxiter=10
+
+		#computational flags: 
+		self.rigid=1
+		self.elastic=1
+		self.eustatic=1
+
+		#numerical discretization accuracy
+		self.degacc=.01
+		
+		#output default:
+		self.requested_outputs=['default']
+
+		#transitions should be a cell array of vectors: 
+		self.transitions=[]
+
+		#default output
+		self.requested_outputs=['default']
+		return self
+		#}}}
+		
+	def checkconsistency(self,md,solution,analyses):    # {{{
+
+		#Early return
+		if (solution!=SealevelriseAnalysisEnum()):
+			return md
+
+		md = checkfield(md,'fieldname','slr.deltathickness','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1])
+		md = checkfield(md,'fieldname','slr.love_h','NaN',1,'Inf',1)
+		md = checkfield(md,'fieldname','slr.love_k','NaN',1,'Inf',1)
+		md = checkfield(md,'fieldname','slr.reltol','size',[1,1])
+		md = checkfield(md,'fieldname','slr.abstol','size',[1,1])
+		md = checkfield(md,'fieldname','slr.maxiter','size',[1,1],'>=',1)
+		md = checkfield(md,'fieldname','slr.degacc','size',[1,1],'>=',1e-10)
+		md = checkfield(md,'fieldname','slr.requested_outputs','stringrow',1)
+
+		#check that love numbers are provided at the same level of accuracy: 
+		if (size(self.love_h,1) != size(self.love_k,1)),
+			error('slr error message: love numbers should be provided at the same level of accuracy')
+		
+		return md
+	# }}}
+	def defaultoutputs(self,md) # {{{
+		return ['SealevelriseS']
+	# }}}
+	def marshall(self,md,fid) # {{{
+		WriteData(fid,'object',self,'class','sealevelrise','fieldname','deltathickness','format','DoubleMat','mattype',1)
+		WriteData(fid,'object',self,'class','sealevelrise','fieldname','reltol','format','Double')
+		WriteData(fid,'object',self,'class','sealevelrise','fieldname','abstol','format','Double')
+		WriteData(fid,'object',self,'class','sealevelrise','fieldname','maxiter','format','Integer')
+		WriteData(fid,'object',self,'class','sealevelrise','fieldname','love_h','format','DoubleMat','mattype',1)
+		WriteData(fid,'object',self,'class','sealevelrise','fieldname','love_k','format','DoubleMat','mattype',1)
+		WriteData(fid,'object',self,'class','sealevelrise','fieldname','rigid','format','Boolean')
+		WriteData(fid,'object',self,'class','sealevelrise','fieldname','elastic','format','Boolean')
+		WriteData(fid,'object',self,'class','sealevelrise','fieldname','eustatic','format','Boolean')
+		WriteData(fid,'object',self,'class','sealevelrise','fieldname','degacc','format','Double')
+		WriteData(fid,'object',self,'class','sealevelrise','fieldname','degacc','format','Double')
+		WriteData(fid,'object',self,'class','sealevelrise','fieldname','transitions','format','MatArray')
+		
+		#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
+		WriteData(fid,'data',outputs,'enum',SealevelriseRequestedOutputsEnum,'format','StringArray')
+
+	# }}}
Index: /issm/trunk-jpl/test/NightlyRun/test2002.py
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test2002.py	(revision 20158)
+++ /issm/trunk-jpl/test/NightlyRun/test2002.py	(revision 20158)
@@ -0,0 +1,94 @@
+#Test Name: EarthSlr
+#Earth Sea Level Rise test. Uses the mesh3dsurface geometry.
+from MatlabFuncs import *
+from model import *
+from EnumDefinitions import *
+from numpy import *
+from parameterize import *
+from solve import *
+
+#mesh earth: 
+md=model() 
+md.mesh=gmshplanet('radius',6.371012*10^3,'resolution',700) #500 km resolution mesh
+
+#parameterize slr solution:
+#slr loading:  {{{
+md.slr.deltathickness=zeros(md.mesh.numberofvertices,1)
+#antarctica
+pos=numpy.nonzero(md.mesh.lat <-80)
+md.slr.deltathickness[pos]=-100
+#greenlnd 
+pos=numpy.nonzero(md.mesh.lat > 70 &  md.mesh.lat < 80 & md.mesh.long>-60 & md.mesh.long<-30)
+md.slr.deltathickness[pos]=-100
+
+#elastic loading from love numbers: 
+love = dlmread('../Data/love_numbers_10k.txt')
+nlov=101
+md.slr.love_h = love(1:nlov,2)  # radial displacement (height) 
+md.slr.love_k = love(1:nlov,4)  # gravitational potential (phi) 
+
+#}}}
+#mask:  {{{
+md.mask=maskpsl() # use maskpsl class (instead of mask) to store the ocean function as a ocean_levelset 
+mask=gmtmask(md.mesh.lat,md.mesh.long) 
+
+icemask=ones(md.mesh.numberofvertices,1)
+pos=find(mask==0)
+icemask(pos)=-1
+pos=find(sum(mask(md.mesh.elements),2)<3)
+icemask(md.mesh.elements(pos))=-1
+
+md.mask.ice_levelset=icemask
+md.mask.ocean_levelset=zeros(md.mesh.numberofvertices,1)
+pos=numpy.nonzero(md.mask.ice_levelset==1)
+md.mask.ocean_levelset(pos)=1
+
+#make sure that the ice level set is all inclusive:
+md.mask.land_levelset=zeros(md.mesh.numberofvertices,1)
+md.mask.groundedice_levelset=-ones(md.mesh.numberofvertices,1) 
+# }}}
+#geometry:  {{{
+di=md.materials.rho_ice/md.materials.rho_water
+md.geometry.thickness=ones(md.mesh.numberofvertices,1)
+md.geometry.surface=(1-di)*zeros(md.mesh.numberofvertices,1)
+md.geometry.base=md.geometry.surface-md.geometry.thickness
+md.geometry.bed=md.geometry.base
+# }}}
+#materials:  {{{
+md.initialization.temperature=273.25*ones(md.mesh.numberofvertices,1)
+md.materials.rheology_B=paterson(md.initialization.temperature)
+md.materials.rheology_n=3*ones(md.mesh.numberofelements,1)
+# }}}
+#Miscellaneous: {{{
+md.miscellaneous.name='slr'
+# }}}
+#Solution parameters:{{{
+md.slr.reltol=NaN
+md.slr.abstol=1e-3
+#}}}
+
+#eustatic run: 
+md.slr.eustatic=1
+md.slr.rigid=0
+md.slr.elastic=0
+md=solve(md,SealevelriseSolutionEnum())
+Seustatic=md.results.SealevelriseSolution.SealevelriseS
+
+#eustatic + rigid run: 
+md.slr.eustatic=1
+md.slr.rigid=1
+md.slr.elastic=0
+md=solve(md,SealevelriseSolutionEnum())
+Srigid=md.results.SealevelriseSolution.SealevelriseS
+
+#eustatic + rigid + elastic run: 
+md.slr.eustatic=1
+md.slr.rigid=1
+md.slr.elastic=1
+md=solve(md,SealevelriseSolutionEnum())
+Selastic=md.results.SealevelriseSolution.SealevelriseS
+
+#Fields and tolerances to track changes
+field_names     =['Eustatic','Rigid','Elastic']
+field_tolerances=[1e-13,1e-13,1e-13]
+field_values=[Seustatic,Srigid,Selastic]