Index: /issm/trunk-jpl/src/m/classes/geometry.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/geometry.m	(revision 26298)
+++ /issm/trunk-jpl/src/m/classes/geometry.m	(revision 26299)
@@ -33,11 +33,4 @@
 	end
 	methods
-		function self = extrude(self,md) % {{{
-			self.surface=project3d(md,'vector',self.surface,'type','node');
-			self.thickness=project3d(md,'vector',self.thickness,'type','node');
-			self.hydrostatic_ratio=project3d(md,'vector',self.hydrostatic_ratio,'type','node');
-			self.base=project3d(md,'vector',self.base,'type','node');
-			self.bed=project3d(md,'vector',self.bed,'type','node');
-		end % }}}
 		function self = geometry(varargin) % {{{
 			switch nargin
@@ -99,4 +92,11 @@
 			WriteData(fid,prefix,'object',self,'fieldname','hydrostatic_ratio','format','DoubleMat','mattype',1);
 		end % }}}
+		function self = extrude(self,md) % {{{
+			self.surface=project3d(md,'vector',self.surface,'type','node');
+			self.thickness=project3d(md,'vector',self.thickness,'type','node');
+			self.hydrostatic_ratio=project3d(md,'vector',self.hydrostatic_ratio,'type','node');
+			self.base=project3d(md,'vector',self.base,'type','node');
+			self.bed=project3d(md,'vector',self.bed,'type','node');
+		end % }}}
 		function savemodeljs(self,fid,modelname) % {{{
 		
Index: /issm/trunk-jpl/src/m/classes/geometry.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/geometry.py	(revision 26298)
+++ /issm/trunk-jpl/src/m/classes/geometry.py	(revision 26299)
@@ -21,6 +21,5 @@
         self.hydrostatic_ratio = np.nan
 
-        nargs = len(args)
-        if nargs == 0:
+        if len(args) == 0:
             self.setdefaultparameters()
         else:
@@ -37,13 +36,4 @@
     #}}}
 
-    def extrude(self, md): #{{{
-        self.surface = project3d(md, 'vector', self.surface, 'type', 'node')
-        self.thickness = project3d(md, 'vector', self.thickness, 'type', 'node')
-        self.hydrostatic_ratio = project3d(md, 'vector', self.hydrostatic_ratio, 'type', 'node')
-        self.base = project3d(md, 'vector', self.base, 'type', 'node')
-        self.bed = project3d(md, 'vector', self.bed, 'type', 'node')
-        return self
-    #}}}
-
     def setdefaultparameters(self): #{{{
         return self
@@ -56,13 +46,13 @@
             md = checkfield(md, 'fieldname', 'geometry.surface', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
             md = checkfield(md, 'fieldname', 'geometry.base', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
-            md = checkfield(md, 'fieldname', 'geometry.thickness', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices], '>', 0, 'timeseries', 1)
-            if any(abs(self.thickness - self.surface + self.base) > 10**-9):
-                md.checkmessage("equality thickness = surface-base violated")
+            md = checkfield(md, 'fieldname', 'geometry.thickness', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices], '>=', 0)
+            if any(abs(self.thickness - self.surface + self.base) > 1e-9):
+                md.checkmessage('equality thickness = surface-base violated')
             if solution == 'TransientSolution' and md.transient.isgroundingline:
                 md = checkfield(md, 'fieldname', 'geometry.bed', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
-                if np.any(self.bed - self.base > 10**-12):
+                if np.any(self.bed - self.base > 1e-12):
                     md.checkmessage('base < bed on one or more vertex')
                 pos = np.where(md.mask.ocean_levelset > 0)
-                if np.any(np.abs(self.bed[pos] - self.base[pos]) > 10**-9):
+                if np.any(np.abs(self.bed[pos] - self.base[pos]) > 1e-9):
                     md.checkmessage('equality base = bed on grounded ice violated')
                 md = checkfield(md, 'fieldname', 'geometry.bed', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
@@ -71,7 +61,8 @@
 
     def marshall(self, prefix, md, fid): #{{{
-        if (len(self.thickness) == md.mesh.numberofvertices) or (len(self.thickness) == md.mesh.numberofvertices + 1):
+        length_thickness = len(self.thickness)
+        if (length_thickness == md.mesh.numberofvertices) or (length_thickness == md.mesh.numberofvertices + 1):
             WriteData(fid, prefix, 'object', self, 'fieldname', 'thickness', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts)
-        elif (len(self.thickness) == md.mesh.numberofelements) or (len(self.thickness) == md.mesh.numberofelements + 1):
+        elif (length_thickness == md.mesh.numberofelements) or (length_thickness == md.mesh.numberofelements + 1):
             WriteData(fid, prefix, 'object', self, 'fieldname', 'thickness', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts)
         else:
@@ -83,2 +74,11 @@
         WriteData(fid, prefix, 'object', self, 'fieldname', 'hydrostatic_ratio', 'format', 'DoubleMat', 'mattype', 1)
     # }}}
+
+    def extrude(self, md): #{{{
+        self.surface = project3d(md, 'vector', self.surface, 'type', 'node')
+        self.thickness = project3d(md, 'vector', self.thickness, 'type', 'node')
+        self.hydrostatic_ratio = project3d(md, 'vector', self.hydrostatic_ratio, 'type', 'node')
+        self.base = project3d(md, 'vector', self.base, 'type', 'node')
+        self.bed = project3d(md, 'vector', self.bed, 'type', 'node')
+        return self
+    #}}}
Index: /issm/trunk-jpl/src/m/classes/materials.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/materials.m	(revision 26298)
+++ /issm/trunk-jpl/src/m/classes/materials.m	(revision 26299)
@@ -96,5 +96,5 @@
 
 					%ice latent heat of fusion L (J/kg)
-					self.latentheat=3.34*10^5;
+					self.latentheat=3.34*1e5;
 
 					%ice thermal conductivity (W/m/K)
@@ -111,5 +111,5 @@
 
 					%rate of change of melting point with pressure (K/Pa)
-					self.beta=9.8*10^-8;
+					self.beta=9.8*1e-8;
 
 					%mixed layer (ice-water interface) heat capacity (J/kg/K)
@@ -117,5 +117,5 @@
 
 					%thermal exchange velocity (ice-water interface) (m/s)
-					self.thermal_exchange_velocity=1.00*10^-4;
+					self.thermal_exchange_velocity=1.00*1e-4;
 
 					%Rheology law: what is the temperature dependence of B with T
@@ -124,5 +124,5 @@
 
 					%Rheology fields default: 
-					self.rheology_B   = 1*1e8; 
+					self.rheology_B   = 1 * 1e8;
 					self.rheology_n   = 3;
 
@@ -136,5 +136,5 @@
 
 					self.viscosity=[1e21;1e40]; %mantle and lithosphere viscosity (respectively) [Pa.s]
-					self.lame_mu=[1.45*1e11;6.7*10^10];  % (Pa) %lithosphere and mantle shear modulus (respectively) [Pa]
+					self.lame_mu=[1.45*1e11;6.7*1e10];  % (Pa) %lithosphere and mantle shear modulus (respectively) [Pa]
 					self.lame_lambda=self.lame_mu;  % (Pa) %mantle and lithosphere lamba parameter (respectively) [Pa]
 					self.burgers_viscosity=[NaN;NaN];
@@ -193,5 +193,5 @@
 				case 'litho'
 					disp(sprintf('   \nLitho:'));
-					fielddisplay(self,'numlayers','number of layers (default 2)');
+					fielddisplay(self,'numlayers','number of layers (default: 2)');
 					fielddisplay(self,'radius','array describing the radius for each interface (numlayers+1) [m]');
 					fielddisplay(self,'viscosity','array describing each layer''s viscosity (numlayers) [Pa.s]');
@@ -207,5 +207,5 @@
 
 
-					fielddisplay(self,'rheologymodel','array describing whether we adopt a MaxWell (0), Burgers (1) or EBM (2) rheology (default 0)');
+					fielddisplay(self,'rheologymodel','array describing whether we adopt a Maxwell (0), Burgers (1) or EBM (2) rheology (default: 0)');
 					fielddisplay(self,'density','array describing each layer''s density (numlayers) [kg/m^3]');
 					fielddisplay(self,'issolid','array describing whether the layer is solid or liquid (default 1) (numlayers)');
@@ -261,9 +261,9 @@
 					end
 					if md.materials.issolid(1)==0 | md.materials.lame_mu(1)==0
-							error('First layer must be solid (issolid(1) > 0 AND lame_mu(1) > 0). Add a weak inner core if necessary.');
+						error('First layer must be solid (issolid(1) > 0 AND lame_mu(1) > 0). Add a weak inner core if necessary.');
 					end
 					ind=find(md.materials.issolid==0);
 					if sum(ismember(diff(ind),1)>=1) %if there are at least two consecutive indices that contain issolid = 0
-							error(['Fluid layers detected at layers #', num2str(ind'), ', but having 2 or more adjacent fluid layers is not supported yet. Consider merging them.'])
+						error(['Fluid layers detected at layers #', num2str(ind'), ', but having 2 or more adjacent fluid layers is not supported yet. Consider merging them.'])
 					end
 
@@ -323,5 +323,5 @@
 					earth_density=0;
 					for i=1:self.numlayers,
-						earth_density=earth_density+ (self.radius(i+1)^3-self.radius(i)^3)*self.density(i);
+						earth_density=earth_density + (self.radius(i+1)^3-self.radius(i)^3)*self.density(i);
 					end
 					earth_density=earth_density/self.radius(self.numlayers+1)^3;
@@ -565,5 +565,5 @@
 	for i=1:length(strnat),
 		switch strnat{i},
-			case 'damageice'
+		case 'damageice'
 			intnat(i)=1;
 		case 'estar'
Index: /issm/trunk-jpl/src/m/classes/materials.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/materials.py	(revision 26298)
+++ /issm/trunk-jpl/src/m/classes/materials.py	(revision 26299)
@@ -16,5 +16,5 @@
     def __init__(self, *args): #{{{
         self.nature = []
-        if not len(args):
+        if len(args) == 0:
             self.nature = ['ice']
         else:
@@ -25,173 +25,188 @@
                 raise RuntimeError("materials constructor error message: nature of the material not supported yet! ('ice' or 'litho' or 'hydro')")
 
-        #start filling in the dynamic fields:
-        for i in range(len(self.nature)):
-            nat = self.nature[i]
-            if nat == 'ice':
-                setattr(self, 'rho_ice', 0)
-                setattr(self, 'rho_water', 0)
-                setattr(self, 'rho_freshwater', 0)
-                setattr(self, 'mu_water', 0)
-                setattr(self, 'heatcapacity', 0)
-                setattr(self, 'latentheat', 0)
-                setattr(self, 'thermalconductivity', 0)
-                setattr(self, 'temperateiceconductivity', 0)
-                setattr(self, 'effectiveconductivity_averaging', 0)
-                setattr(self, 'meltingpoint', 0)
-                setattr(self, 'beta', 0)
-                setattr(self, 'mixed_layer_capacity', 0)
-                setattr(self, 'thermal_exchange_velocity', 0)
-                setattr(self, 'rheology_B', 0)
-                setattr(self, 'rheology_n', 0)
-                setattr(self, 'rheology_law', 0)
-            elif nat == 'litho':
-                setattr(self, 'numlayers', 0)
-                setattr(self, 'radius', 0)
-                setattr(self, 'viscosity', 0)
-                setattr(self, 'lame_lambda', 0)
-                setattr(self, 'lame_mu', 0)
-                setattr(self, 'burgers_viscosity', 0)
-                setattr(self, 'burgers_mu', 0)
-                setattr(self, 'isburgers', 0)
-                setattr(self, 'density', 0)
-                setattr(self, 'issolid', 0)
-            elif nat == 'hydro':
-                setattr(self, 'rho_ice', 0)
-                setattr(self, 'rho_water', 0)
-                setattr(self, 'rho_freshwater', 0)
+        for i in range(len(self.nature)):
+            nat = self.nature[i]
+            if nat == 'ice':
+                self.rho_ice = 0
+                self.rho_water = 0
+                self.rho_freshwater = 0
+                self.mu_water = 0
+                self.heatcapacity = 0
+                self.latentheat = 0
+                self.thermalconductivity = 0
+                self.temperateiceconductivity = 0
+                self.effectiveconductivity_averaging = 0
+                self.meltingpoint = 0
+                self.beta = 0
+                self.mixed_layer_capacity = 0
+                self.thermal_exchange_velocity = 0
+                self.rheology_B = 0
+                self.rheology_n = 0
+                self.rheology_law = 0
+            elif nat == 'litho':
+                self.numlayers = 0
+                self.radius = 0
+                self.viscosity = 0
+                self.lame_lambda = 0
+                self.lame_mu = 0
+                self.burgers_viscosity = 0
+                self.burgers_mu = 0
+                self.ebm_alpha = 0
+                self.ebm_delta = 0
+                self.ebm_taul = 0
+                self.ebm_tauh = 0
+                self.rheologymodel = 0
+                self.density = 0
+                self.issolid = 0
+            elif nat == 'hydro':
+                self.rho_ice = 0
+                self.rho_water = 0
+                self.rho_freshwater = 0
             else:
                 raise RuntimeError("materials constructor error message: nature of the material not supported yet! ('ice' or 'litho' or 'hydro')")
-        setattr(self, 'earth_density', 0)
-
-        #set default parameters:
+        self.earth_density = 0
+
+        # Set default parameters
         self.setdefaultparameters()
     #}}}
 
+    def __repr__(self): #{{{
+        s = '   Materials:\n'
+        for i in range(len(self.nature)):
+            nat = self.nature[i]
+            if nat == 'ice':
+                s += 'Ice:\n'
+                s += '{}\n'.format(fielddisplay(self, 'rho_ice', 'ice density [kg/m^3]'))
+                s += '{}\n'.format(fielddisplay(self, 'rho_water', 'ocean water density [kg/m^3]'))
+                s += '{}\n'.format(fielddisplay(self, 'rho_freshwater', 'fresh water density [kg/m^3]'))
+                s += '{}\n'.format(fielddisplay(self, 'mu_water', 'water viscosity [N s/m^2]'))
+                s += '{}\n'.format(fielddisplay(self, 'heatcapacity', 'heat capacity [J/kg/K]'))
+                s += '{}\n'.format(fielddisplay(self, 'thermalconductivity', 'ice thermal conductivity [W/m/K]'))
+                s += '{}\n'.format(fielddisplay(self, 'temperateiceconductivity', 'temperate ice thermal conductivity [W/m/K]'))
+                s += '{}\n'.format(fielddisplay(self, 'meltingpoint', 'melting point of ice at 1atm in K'))
+                s += '{}\n'.format(fielddisplay(self, 'latentheat', 'latent heat of fusion [J/m^3]'))
+                s += '{}\n'.format(fielddisplay(self, 'beta', 'rate of change of melting point with pressure [K/Pa]'))
+                s += '{}\n'.format(fielddisplay(self, 'mixed_layer_capacity', 'mixed layer capacity [W/kg/K]'))
+                s += '{}\n'.format(fielddisplay(self, 'thermal_exchange_velocity', 'thermal exchange velocity [m/s]'))
+                s += '{}\n'.format(fielddisplay(self, 'rheology_B', 'flow law parameter [Pa s^(1/n)]'))
+                s += '{}\n'.format(fielddisplay(self, 'rheology_n', 'Glen\'s flow law exponent'))
+                s += '{}\n'.format(fielddisplay(self, 'rheology_law', 'law for the temperature dependance of the rheology: \'None\', \'BuddJacka\', \'Cuffey\', \'CuffeyTemperate\', \'Paterson\', \'Arrhenius\', \'LliboutryDuval\', \'NyeCO2\', or \'NyeH2O\''))
+            elif nat == 'litho':
+                s += 'Litho:\n'
+                s += '{}\n'.format(fielddisplay(self, 'numlayers', 'number of layers (default: 2)'))
+                s += '{}\n'.format(fielddisplay(self, 'radius', 'array describing the radius for each interface (numlayers + 1) [m]'))
+                s += '{}\n'.format(fielddisplay(self, 'viscosity', 'array describing each layer\'s viscosity (numlayers) [Pa.s]'))
+                s += '{}\n'.format(fielddisplay(self, 'lame_lambda', 'array describing the lame lambda parameter (numlayers) [Pa]'))
+                s += '{}\n'.format(fielddisplay(self, 'lame_mu', 'array describing the shear modulus for each layers (numlayers) [Pa]'))
+                s += '{}\n'.format(fielddisplay(self, 'burgers_viscosity', 'array describing each layer\'s transient viscosity, only for Burgers rheologies  (numlayers) [Pa.s]'))
+                s += '{}\n'.format(fielddisplay(self, 'burgers_mu', 'array describing each layer\'s transient shear modulus, only for Burgers rheologies  (numlayers) [Pa]'))
+
+                s += '{}\n'.format(fielddisplay(self, 'ebm_alpha', 'array describing each layer\'s exponent parameter controlling the shape of shear modulus curve between taul and tauh, only for EBM rheology (numlayers)'))
+                s += '{}\n'.format(fielddisplay(self, 'ebm_delta', 'array describing each layer\'s amplitude of the transient relaxation (ratio between elastic rigity to pre-maxwell relaxation rigity), only for EBM rheology (numlayers)'))
+                s += '{}\n'.format(fielddisplay(self, 'ebm_taul', 'array describing each layer\'s starting period for transient relaxation, only for EBM rheology  (numlayers) [s]'))
+                s += '{}\n'.format(fielddisplay(self, 'ebm_tauh', 'array describing each layer''s array describing each layer\'s end period for transient relaxation, only for Burgers rheology  (numlayers) [s]'))
+
+                s += '{}\n'.format(fielddisplay(self, 'rheologymodel', 'array describing whether we adopt a Maxwell (0), Burgers (1) or EBM (2) rheology (default: 0)'))
+                s += '{}\n'.format(fielddisplay(self, 'density', 'array describing each layer\'s density (numlayers) [kg/m^3]'))
+                s += '{}\n'.format(fielddisplay(self, 'issolid', 'array describing whether the layer is solid or liquid (default 1) (numlayers)'))
+            elif nat == 'hydro':
+                s += 'Hydro:\n'
+                s += '{}\n'.format(fielddisplay(self, 'rho_ice', 'ice density [kg/m^3]'))
+                s += '{}\n'.format(fielddisplay(self, 'rho_water', 'ocean water density [kg/m^3]'))
+                s += '{}\n'.format(fielddisplay(self, 'earth_density', 'mantle density [kg/m^3]'))
+                s += '{}\n'.format(fielddisplay(self, 'rho_freshwater', 'fresh water density [kg/m^3]'))
+
+            else:
+                raise RuntimeError('materials constructor error message: nature of the material not supported yet! (\'ice\' or \'litho\' or \'hydro\')')
+        return s
+    #}}}
+
     def setdefaultparameters(self): #{{{
         for i in range(len(self.nature)):
             nat = self.nature[i]
             if nat == 'ice':
-                #ice density (kg/m^3)
+                # Ice density (kg/m^3)
                 self.rho_ice = 917.
 
-                #ocean water density (kg/m^3)
+                # Ocean water density (kg/m^3)
                 self.rho_water = 1023.
 
-                #fresh water density (kg/m^3)
+                # Fresh water density (kg/m^3)
                 self.rho_freshwater = 1000.
 
-                #water viscosity (N.s/m^2)
+                # Water viscosity (N.s/m^2)
                 self.mu_water = 0.001787
 
-                #ice heat capacity cp (J/kg/K)
+                # Ice heat capacity cp (J/kg/K)
                 self.heatcapacity = 2093.
 
-                #ice latent heat of fusion L (J / kg)
-                self.latentheat = 3.34e5
-
-                #ice thermal conductivity (W/m/K)
+                # Ice latent heat of fusion L (J/kg)
+                self.latentheat = 3.34 * 1e5
+
+                # Ice thermal conductivity (W/m/K)
                 self.thermalconductivity = 2.4
 
-                #wet ice thermal conductivity (W/m/K)
+                # Wet ice thermal conductivity (W/m/K)
                 self.temperateiceconductivity = 0.24
 
-                #computation of effective conductivity
+                # Computation of effective conductivity
                 self.effectiveconductivity_averaging = 1
 
-                #the melting point of ice at 1 atmosphere of pressure in K
+                # The melting point of ice at 1 atmosphere of pressure in K
                 self.meltingpoint = 273.15
 
-                #rate of change of melting point with pressure (K/Pa)
-                self.beta = 9.8e-8
-
-                #mixed layer (ice-water interface) heat capacity (J/kg/K)
+                # Rate of change of melting point with pressure (K/Pa)
+                self.beta = 9.8 * 1e-8
+
+                # Mixed layer (ice-water interface) heat capacity (J/kg/K)
                 self.mixed_layer_capacity = 3974.
 
-                #thermal exchange velocity (ice-water interface) (m/s)
-                self.thermal_exchange_velocity = 1.00e-4
-
-                #Rheology law: what is the temperature dependence of B with T
-                #available: none, paterson and arrhenius
+                # Thermal exchange velocity (ice-water interface) (m/s)
+                self.thermal_exchange_velocity = 1.00 * 1e-4
+
+                # Rheology law: what is the temperature dependence of B with T
+                # available: none, paterson and arrhenius
                 self.rheology_law = 'Paterson'
 
-                #Rheology fields default
-                self.rheology_B = 1e8
+                # Rheology fields default
+                self.rheology_B = 1 * 1e8
                 self.rheology_n = 3
             elif nat == 'litho':
-                #we default to a configuration that enables running GIA solutions using giacaron and/or giaivins.
+                # We default to a configuration that enables running GIA 
+                # solutions using giacaron and/or giaivins
                 self.numlayers = 2
 
-                #center of the earth (approximation, must not be 0), then the lab (lithosphere / asthenosphere boundary) then the surface
-                #(with 1d3 to avoid numerical singularities)
-                self.radius = [1e3, 6278e3, 6378e3]
-
-                self.viscosity = [1e21, 1e40]  #mantle and lithosphere viscosity (respectively) [Pa.s]
-                self.lame_mu = [1.45e11, 6.7e10]  # (Pa)  #lithosphere and mantle shear modulus (respectively) [Pa]
-                self.lame_lambda = self.lame_mu  # (Pa)  #mantle and lithosphere lamba parameter (respectively) [Pa]
+                # Center of the earth (approximation, must not be 0), then the 
+                # lab (lithosphere/asthenosphere boundary) then the surface
+                # (with 1d3 to avoid numerical singularities)
+                self.radius = [1e3, 6278 * 1e3, 6378 * 1e3]
+
+                self.viscosity = [1e21, 1e40] # Mantle and lithosphere viscosity (respectively) [Pa.s]
+                self.lame_mu = [1.45 * 1e11, 6.7 * 1e10] # (Pa) # Lithosphere and mantle shear modulus (respectively) [Pa]
+                self.lame_lambda = self.lame_mu # (Pa) # Mantle and lithosphere lamba parameter (respectively) [Pa]
                 self.burgers_viscosity = [np.nan, np.nan]
                 self.burgers_mu = [np.nan, np.nan]
-                self.isburgers = [0, 0]
-                self.density = [5.51e3, 5.50e3]  # (Pa)  #mantle and lithosphere density [kg/m^3]
-                self.issolid = [1, 1]  # is layer solid or liquid.
-            elif nat == 'hydro':
-                #ice density (kg/m^3)
+
+                self.ebm_alpha = [np.nan, np.nan]
+                self.ebm_delta = [np.nan, np.nan]
+                self.ebm_taul = [np.nan, np.nan]
+                self.ebm_tauh = [np.nan, np.nan]
+                self.rheologymodel = [0, 0]
+                self.density = [5.51e3, 5.50e3]  # (Pa) # Mantle and lithosphere density [kg/m^3]
+                self.issolid = [1, 1]  # Is layer solid or liquid?
+            elif nat == 'hydro':
+                # Ice density (kg/m^3)
                 self.rho_ice = 917.
 
-                #ocean water density (kg/m^3)
+                # Ocean water density (kg/m^3)
                 self.rho_water = 1023.
 
-                #fresh water density (kg/m^3)
+                # Fresh water density (kg/m^3)
                 self.rho_freshwater = 1000.
             else:
                 raise RuntimeError("materials setdefaultparameters error message: nature of the material not supported yet! ('ice' or 'litho' or 'hydro')")
 
-            #average density of the Earth (kg/m^3)
+            # Average density of the Earth (kg/m^3)
             self.earth_density = 5512
-    #}}}
-
-    def __repr__(self): #{{{
-        string = "   Materials:"
-        for i in range(len(self.nature)):
-            nat = self.nature[i]
-            if nat == 'ice':
-                string = "%s\n%s" % (string, 'Ice:')
-                string = "%s\n%s" % (string, fielddisplay(self, "rho_ice", "ice density [kg/m^3]"))
-                string = "%s\n%s" % (string, fielddisplay(self, "rho_water", "ocean water density [kg/m^3]"))
-                string = "%s\n%s" % (string, fielddisplay(self, "rho_freshwater", "fresh water density [kg/m^3]"))
-                string = "%s\n%s" % (string, fielddisplay(self, "mu_water", "water viscosity [N s/m^2]"))
-                string = "%s\n%s" % (string, fielddisplay(self, "heatcapacity", "heat capacity [J/kg/K]"))
-                string = "%s\n%s" % (string, fielddisplay(self, "thermalconductivity", "ice thermal conductivity [W/m/K]"))
-                string = "%s\n%s" % (string, fielddisplay(self, "temperateiceconductivity", "temperate ice thermal conductivity [W/m/K]"))
-                string = "%s\n%s" % (string, fielddisplay(self, "meltingpoint", "melting point of ice at 1atm in K"))
-                string = "%s\n%s" % (string, fielddisplay(self, "latentheat", "latent heat of fusion [J/m^3]"))
-                string = "%s\n%s" % (string, fielddisplay(self, "beta", "rate of change of melting point with pressure [K/Pa]"))
-                string = "%s\n%s" % (string, fielddisplay(self, "mixed_layer_capacity", "mixed layer capacity [W/kg/K]"))
-                string = "%s\n%s" % (string, fielddisplay(self, "thermal_exchange_velocity", "thermal exchange velocity [m/s]"))
-                string = "%s\n%s" % (string, fielddisplay(self, "rheology_B", "flow law parameter [Pa s^(1/n)]"))
-                string = "%s\n%s" % (string, fielddisplay(self, "rheology_n", "Glen's flow law exponent"))
-                string = "%s\n%s" % (string, fielddisplay(self, "rheology_law", "law for the temperature dependance of the rheology: 'None', 'BuddJacka', 'Cuffey', 'CuffeyTemperate', 'Paterson', 'Arrhenius', 'LliboutryDuval', 'NyeCO2', or 'NyeH2O'"))
-            elif nat == 'litho':
-                string = "%s\n%s" % (string, 'Litho:')
-                string = "%s\n%s" % (string, fielddisplay(self, 'numlayers', 'number of layers (default 2)'))
-                string = "%s\n%s" % (string, fielddisplay(self, 'radius', 'array describing the radius for each interface (numlayers + 1) [m]'))
-                string = "%s\n%s" % (string, fielddisplay(self, 'viscosity', 'array describing each layer''s viscosity (numlayers) [Pa.s]'))
-                string = "%s\n%s" % (string, fielddisplay(self, 'lame_lambda', 'array describing the lame lambda parameter (numlayers) [Pa]'))
-                string = "%s\n%s" % (string, fielddisplay(self, 'lame_mu', 'array describing the shear modulus for each layers (numlayers) [Pa]'))
-                string = "%s\n%s" % (string, fielddisplay(self, 'burgers_viscosity', 'array describing each layer''s transient viscosity, only for Burgers rheologies  (numlayers) [Pa.s]'))
-                string = "%s\n%s" % (string, fielddisplay(self, 'burgers_mu', 'array describing each layer''s transient shear modulus, only for Burgers rheologies  (numlayers) [Pa]'))
-                string = "%s\n%s" % (string, fielddisplay(self, 'isburgers', 'array describing whether we adopt a MaxWell (0) or Burgers (1) rheology (default 0)'))
-                string = "%s\n%s" % (string, fielddisplay(self, 'density', 'array describing each layer''s density (numlayers) [kg/m^3]'))
-                string = "%s\n%s" % (string, fielddisplay(self, 'issolid', 'array describing whether the layer is solid or liquid (default 1) (numlayers)'))
-            elif nat == 'hydro':
-                string = "%s\n%s" % (string, 'Hydro:')
-                string = "%s\n%s" % (string, fielddisplay(self, "rho_ice", "ice density [kg/m^3]"))
-                string = "%s\n%s" % (string, fielddisplay(self, "rho_water", "ocean water density [kg/m^3]"))
-                string = "%s\n%s" % (string, fielddisplay(self, "earth_density", "mantle density [kg/m^3]"))
-                string = "%s\n%s" % (string, fielddisplay(self, "rho_freshwater", "fresh water density [kg/m^3]"))
-
-            else:
-                raise RuntimeError("materials constructor error message: nature of the material not supported yet! ('ice' or 'litho' or 'hydro')")
-
-        return string
     #}}}
 
@@ -218,18 +233,24 @@
                 md = checkfield(md, 'fieldname', 'materials.density', 'NaN', 1, 'Inf', 1, 'size', [md.materials.numlayers, 1], '>', 0)
                 md = checkfield(md, 'fieldname', 'materials.viscosity', 'NaN', 1, 'Inf', 1, 'size', [md.materials.numlayers, 1], '>=', 0)
-                md = checkfield(md, 'fieldname', 'materials.isburgers', 'NaN', 1, 'Inf', 1, 'size', [md.materials.numlayers, 1], '>=', 0, '<=', 1)
+                md = checkfield(md, 'fieldname', 'materials.rheologymodel', 'NaN', 1, 'Inf', 1, 'size', [md.materials.numlayers, 1], '>=', 0, '<=', 2)
                 md = checkfield(md, 'fieldname', 'materials.burgers_viscosity', 'Inf', 1, 'size', [md.materials.numlayers, 1], '>=', 0)
                 md = checkfield(md, 'fieldname', 'materials.burgers_mu', 'Inf', 1, 'size', [md.materials.numlayers, 1], '>=', 0)
+                md = checkfield(md, 'fieldname', 'materials.ebm_alpha', 'Inf', 1, 'size', [md.materials.numlayers, 1], '>=', 0)
+                md = checkfield(md, 'fieldname', 'materials.ebm_delta', 'Inf', 1, 'size', [md.materials.numlayers, 1], '>=', 0)
+                md = checkfield(md, 'fieldname', 'materials.ebm_taul', 'Inf', 1, 'size', [md.materials.numlayers, 1], '>=', 0)
+                md = checkfield(md, 'fieldname', 'materials.ebm_tauh', 'Inf', 1, 'size', [md.materials.numlayers, 1], '>=', 0)
 
                 for i in range(md.materials.numlayers):
-                    if md.materials.isburgers[i] and (np.isnan(md.materials.burgers_viscosity[i] or np.isnan(md.materials.burgers_mu[i]))):
-                        raise RuntimeError("materials checkconsistency error message: Litho burgers_viscosity or burgers_mu has NaN values, inconsistent with isburgers choice")
-
-                    if md.materials.issolid[0] == 0 or md.materials.lame_mu[0] == 0:
-                        raise RuntimeError('First layer must be solid (issolid(1) > 0 AND lame_mu(1) > 0). Add a weak inner core if necessary.')
-
-                    for i in range(md.materials.numlayers - 1):
-                        if (not md.materials.issolid[i]) and (not md.materials.issolid[i + 1]):  #if there are at least two consecutive indices that contain issolid = 0
-                            raise RuntimeError("%s%i%s" % ('2 or more adjacent fluid layers detected starting at layer ', i, '. This is not supported yet. Consider merging them.'))
+                    if md.materials.rheologymodel[i] == 1 and (np.isnan(md.materials.burgers_viscosity[i] or np.isnan(md.materials.burgers_mu[i]))):
+                        raise RuntimeError('materials checkconsistency error message: Litho burgers_viscosity or burgers_mu has NaN values, inconsistent with rheologymodel choice')
+
+                    if md.materials.rheologymodel[i] == 2 and (np.isnan(md.materials.ebm_alpha[i]) or np.isnan(md.materials.ebm_delta[i]) or np.isnan(md.materials.ebm_taul[i]) or np.isnan(md.materials.ebm_tauh[i])):
+                        raise RuntimeError('materials checkconsistency error message: Litho ebm_alpha, ebm_delta, ebm_taul or ebm_tauh has NaN values, inconsistent with rheologymodel choice')
+                if md.materials.issolid[0] == 0 or md.materials.lame_mu[0] == 0:
+                    raise RuntimeError('First layer must be solid (issolid[0] > 0 AND lame_mu[0] > 0). Add a weak inner core if necessary.')
+                ind = np.where(md.materials.issolid == 0)[0]
+                if np.sum(np.in1d(np.diff(ind),1) >= 1): # If there are at least two consecutive indices that contain issolid = 0
+                    raise RuntimeError('Fluid layers detected at layers #{} but having 2 or more adjacent fluid layers is not supported yet. Consider merging them.'.format(i))
+
             elif nat == 'hydro':
                 md = checkfield(md, 'fieldname', 'materials.rho_ice', '>', 0)
@@ -238,5 +259,5 @@
                 md = checkfield(md, 'fieldname', 'materials.rho_freshwater', '>', 0)
             else:
-                raise RuntimeError("materials checkconsistency error message: nature of the material not supported yet! ('ice' or 'litho' or 'hydro')")
+                raise RuntimeError('materials checkconsistency error message: nature of the material not supported yet! (\'ice\' or \'litho\' or \'hydro\')')
 
         return md
@@ -250,5 +271,4 @@
             nat = self.nature[i]
             if nat == 'ice':
-                WriteData(fid, prefix, 'name', 'md.materials.type', 'data', 3, 'format', 'Integer')
                 WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'rho_ice', 'format', 'Double')
                 WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'rho_water', 'format', 'Double')
@@ -275,12 +295,16 @@
                 WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'density', 'format', 'DoubleMat', 'mattype', 3)
                 WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'viscosity', 'format', 'DoubleMat', 'mattype', 3)
-                WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'isburgers', 'format', 'DoubleMat', 'mattype', 3)
+                WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'rheologymodel', 'format', 'DoubleMat', 'mattype', 3) 
                 WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'burgers_viscosity', 'format', 'DoubleMat', 'mattype', 3)
                 WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'burgers_mu', 'format', 'DoubleMat', 'mattype', 3)
+                WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'ebm_alpha', 'format', 'DoubleMat', 'mattype', 3)
+                WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'ebm_delta', 'format', 'DoubleMat', 'mattype', 3)
+                WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'ebm_taul', 'format', 'DoubleMat', 'mattype', 3)
+                WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'ebm_tauh', 'format', 'DoubleMat', 'mattype', 3)
                 # Compute earth density compatible with our layer density distribution
                 earth_density = 0
-                for i in range(len(self.numlayers)):
-                    earth_density = earth_density + (self.radius[i + 1] ** 3 - self.radius[i] ** 3) * self.density[i]
-                earth_density = earth_density / self.radius[self.numlayers + 1] ** 3
+                for i in range(self.numlayers):
+                    earth_density = earth_density + (pow(self.radius[i + 1], 3) - pow(self.radius[i], 3)) * self.density[i]
+                earth_density = earth_density / pow(self.radius[self.numlayers], 3)
                 self.earth_density = earth_density
             elif nat == 'hydro':
@@ -289,5 +313,5 @@
                 WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'rho_freshwater', 'format', 'Double')
             else:
-                raise RuntimeError("materials constructor error message: nature of the material not supported yet! ('ice' or 'litho' or 'hydro')")
+                raise RuntimeError('materials constructor error message: nature of the material not supported yet! (\'ice\' or \'litho\' or \'hydro\')')
         WriteData(fid, prefix, 'data', self.earth_density, 'name', 'md.materials.earth_density', 'format', 'Double')
     #}}}
@@ -322,5 +346,5 @@
             intnat[i] = 7
         else:
-            raise RuntimeError("materials constructor error message: nature of the material not supported yet! ('ice' or 'litho' or 'hydro')")
+            raise RuntimeError('materials constructor error message: nature of the material not supported yet! (\'ice\' or \'litho\' or \'hydro\')')
 
     return intnat
Index: /issm/trunk-jpl/src/m/classes/solidearth.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/solidearth.m	(revision 26298)
+++ /issm/trunk-jpl/src/m/classes/solidearth.m	(revision 26299)
@@ -83,5 +83,4 @@
 			end
 
-			
 		end % }}}
 		function list=defaultoutputs(self,md) % {{{
@@ -107,8 +106,8 @@
 		end % }}}
 		function marshall(self,prefix,md,fid) % {{{
-			
+
 			WriteData(fid,prefix,'object',self,'fieldname','planetradius','format','Double');
 			WriteData(fid,prefix,'object',self,'fieldname','transitions','format','MatArray');
-		
+
 			if ~isempty(self.partitionice),
 				npartice=max(self.partitionice)+2;
@@ -127,6 +126,4 @@
 			end
 
-
-			
 			WriteData(fid,prefix,'object',self,'fieldname','partitionice','mattype',1,'format','DoubleMat');
 			WriteData(fid,prefix,'data',npartice,'format','Integer','name','md.solidearth.npartice');
@@ -135,5 +132,5 @@
 			WriteData(fid,prefix,'object',self,'fieldname','partitionocean','mattype',1,'format','DoubleMat');
 			WriteData(fid,prefix,'data',npartocean,'format','Integer','name','md.solidearth.npartocean');
-			
+
 			self.settings.marshall(prefix,md,fid);
 			self.lovenumbers.marshall(prefix,md,fid);
@@ -145,8 +142,8 @@
 				WriteData(fid,prefix,'data',0,'format','Integer','name','md.solidearth.isexternal');
 			end
-			
+
 			%process requested outputs
 			outputs = self.requested_outputs;
-			pos  = find(ismember(outputs,'default'));
+			pos = find(ismember(outputs,'default'));
 			if ~isempty(pos),
 				outputs(pos) = [];                         %remove 'default' from outputs
@@ -155,4 +152,6 @@
 			WriteData(fid,prefix,'data',outputs,'name','md.solidearth.requested_outputs','format','StringArray');
 
+		end % }}}
+		function self = extrude(self,md) % {{{
 		end % }}}
 		function savemodeljs(self,fid,modelname) % {{{
@@ -168,6 +167,4 @@
 			writejscellarray(fid,[modelname '.solidearth.partition'],self.partition);
 		end % }}}
-		function self = extrude(self,md) % {{{
-		end % }}}
 	end
 end
Index: /issm/trunk-jpl/src/m/classes/solidearth.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/solidearth.py	(revision 26298)
+++ /issm/trunk-jpl/src/m/classes/solidearth.py	(revision 26299)
@@ -19,9 +19,7 @@
     """
 
-    def __init__(self, *args):  #{{{
-        self.initialsealevel    = np.nan
+    def __init__(self, *args): #{{{
         self.settings           = solidearthsettings()
         self.external           = []
-        self.surfaceload        = surfaceload()
         self.lovenumbers        = lovenumbers()
         self.rotational         = rotational()
@@ -31,7 +29,7 @@
         self.partitionice       = []
         self.partitionhydro     = []
+        self.partitionocean     = []
 
         nargs = len(args)
-
         if nargs == 0:
             self.setdefaultparameters('earth')
@@ -41,7 +39,6 @@
             raise Exception('solidearth constructor error message: zero or one argument only!')
     #}}}
-    def __repr__(self):  # {{{
+    def __repr__(self): #{{{
         s = '   solidearthinputs, forcings and settings:\n'
-        s += '{}\n'.format(fielddisplay(self, 'initialsealevel', 'sea level at the start of computation [m]'))
         s += '{}\n'.format(fielddisplay(self, 'planetradius', 'planet radius [m]'))
         s += '{}\n'.format(fielddisplay(self, 'transitions', 'indices into parts of the mesh that will be icecaps'))
@@ -49,8 +46,8 @@
         s += '{}\n'.format(fielddisplay(self, 'partitionice', 'ice partition vector for barystatic contribution'))
         s += '{}\n'.format(fielddisplay(self, 'partitionhydro', 'hydro partition vector for barystatic contribution'))
+        s += '{}\n'.format(fielddisplay(self, 'partitionocean', 'ocean partition vector for barystatic contribution'))
         if not self.external:
             s += '{}\n'.format(fielddisplay(self, 'external', 'external solution, of the type solidearthsolution'))
         print(self.settings)
-        print(self.surfaceload)
         print(self.lovenumbers)
         print(self.rotational)
@@ -59,6 +56,6 @@
         return s
     #}}}
-    def setdefaultparameters(self, planet):  # {{{
-        # Default output
+    def setdefaultparameters(self, planet): #{{{
+        # Output default
         self.requested_outputs = ['default']
 
@@ -69,6 +66,7 @@
         self.partitionice = []
         self.partitionhydro = []
+        self.partitionocean = []
 
-        # No external solutions by defalt
+        # No external solutions by default
         self.external = []
 
@@ -76,26 +74,24 @@
         self.planetradius = planetradius(planet)
     #}}}
-    def checkconsistency(self, md, solution, analyses):  # {{{
+    def checkconsistency(self, md, solution, analyses): #{{{
         if ('SealevelriseAnalysis' not in analyses) or (solution == 'TransientSolution' and not md.transient.isslc):
             return md
-        md = checkfield(md, 'fieldname', 'solidearth.initialsealevel', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
+
         md = checkfield(md, 'fieldname', 'solidearth.requested_outputs', 'stringrow', 1)
 
         self.settings.checkconsistency(md, solution, analyses)
-        self.surfaceload.checkconsistency(md, solution, analyses)
         self.lovenumbers.checkconsistency(md, solution, analyses)
         self.rotational.checkconsistency(md, solution, analyses)
         if self.external:
-            if not isinstance(self.external,'solidearthsolution'):
+            if not isinstance(self.external, 'solidearthsolution'):
                 raise Exception('solidearth consistency check: external field should be a solidearthsolution')
             end
-            self.external.checkconsistency(md,solution,analyses)
+            self.external.checkconsistency(md, solution, analyses)
         return md
     #}}}
-    def defaultoutputs(self, md):  #{{{
+    def defaultoutputs(self, md): #{{{
         return ['Sealevel']
     #}}}
-    def marshall(self, prefix, md, fid):  #{{{
-        WriteData(fid, prefix, 'object', self, 'fieldname', 'initialsealevel', 'mattype', 1, 'format', 'DoubleMat', 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts)
+    def marshall(self, prefix, md, fid): #{{{
         WriteData(fid, prefix, 'object', self, 'fieldname', 'planetradius', 'format', 'Double')
         WriteData(fid, prefix, 'object', self, 'fieldname', 'transitions', 'format', 'MatArray')
@@ -111,15 +107,24 @@
             nparthydro = 0
 
+        if len(self.partitionocean):
+            npartocean = np.max(self.partitionocean) + 2
+        else:
+            npartocean = 0
+
         WriteData(fid, prefix, 'object', self, 'fieldname', 'partitionice', 'mattype', 1, 'format', 'DoubleMat');
         WriteData(fid, prefix, 'data', npartice, 'format', 'Integer', 'name', 'md.solidearth.npartice');
         WriteData(fid, prefix, 'object', self, 'fieldname', 'partitionhydro', 'mattype', 1, 'format', 'DoubleMat');
         WriteData(fid, prefix, 'data', nparthydro,'format', 'Integer', 'name','md.solidearth.nparthydro');
+        WriteData(fid, prefix, 'object', self, 'fieldname', 'partitionocean', 'mattype', 1, 'format', 'DoubleMat');
+        WriteData(fid, prefix, 'data', npartocean,'format', 'Integer', 'name','md.solidearth.npartocean');
 
         self.settings.marshall(prefix, md, fid)
-        self.surfaceload.marshall(prefix, md, fid)
         self.lovenumbers.marshall(prefix, md, fid)
         self.rotational.marshall(prefix, md, fid)
         if self.external:
+            WriteData(fid, prefix, 'data', 1, 'format', 'Integer', 'name', 'md.solidearth.isexternal')
             self.external.marshall(prefix, md, fid)
+        else:
+            WriteData(fid, prefix, 'data', 0, 'format', 'Integer', 'name', 'md.solidearth.isexternal')
 
         #process requested outputs
@@ -127,10 +132,9 @@
         pos = np.where(np.asarray(outputs) == 'default')[0]
         if len(pos):
-            outputs = np.delete(outputs, pos)  #remove 'default' from outputs
-            outputs = np.append(outputs, self.defaultoutputs(md))  #add defaults
+            outputs = np.delete(outputs, pos)  # remove 'default' from outputs
+            outputs = np.append(outputs, self.defaultoutputs(md)) # add defaults
         WriteData(fid, prefix, 'data', outputs, 'name', 'md.solidearth.requested_outputs', 'format', 'StringArray')
     #}}}
     def extrude(self, md): #{{{
-        self.initialsealevel = project3d(md, 'vector', self.initialsealevel, 'type', 'node')
         return self
     #}}}
Index: /issm/trunk-jpl/src/m/classes/solidearthsettings.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/solidearthsettings.m	(revision 26298)
+++ /issm/trunk-jpl/src/m/classes/solidearthsettings.m	(revision 26299)
@@ -5,5 +5,5 @@
 
 classdef solidearthsettings
-	properties (SetAccess=public) 
+	properties (SetAccess=public)
 		reltol                 = 0;
 		abstol                 = 0;
@@ -13,5 +13,5 @@
 		viscous                = 1;
 		rotation               = 1;
-		grdocean	       = 1;
+		grdocean               = 1;
 		ocean_area_scaling     = 0;
 		runfrequency           = 1; %how many time steps we skip before we run grd_core
@@ -20,7 +20,7 @@
 		compute_bp_grd         = 0; %will GRD patterns for bottom pressures be computed? 
 		degacc                 = 0; %degree increment for resolution of Green tables.
-		timeacc                = 0; %time step accurary required to compute Green tables
+		timeacc                = 0; %time step accuracy required to compute Green tables
 		horiz                  = 0; %compute horizontal deformation
-		grdmodel               = 0; %grd model (0 by default, 1 for (visco)-elastic, 2 for Ivins)
+		grdmodel               = 0; %grd model (0 by default, 1 for (visco-)elastic, 2 for Ivins)
 		cross_section_shape    = 0; %cross section only used when grd model is Ivins
 	end
@@ -81,9 +81,9 @@
 		%how many time steps we skip before we run solidearthsettings solver during transient
 		self.runfrequency=1;
-		
-		%horizontal displacement?  (not by default)
+
+		%horizontal displacement? (not on by default)
 		self.horiz=0;
-		
-		%cross section for Ivins model 
+
+		%cross section for Ivins model
 		self.cross_section_shape=1; %square as default (see iedge in GiaDeflectionCorex)
 
@@ -106,4 +106,5 @@
 			md = checkfield(md,'fieldname','solidearth.settings.grdmodel','>=',0,'<=',2);
 			md = checkfield(md,'fieldname','solidearth.settings.cross_section_shape','numel',[1],'values',[1,2]);
+
 			if self.elastic==1 & self.selfattraction==0,
 				error('solidearthsettings checkconsistency error message: need selfattraction on if elastic flag is set');
@@ -137,12 +138,12 @@
 			disp(sprintf('   solidearth settings:'));
 
-			fielddisplay(self,'reltol','sea level change relative convergence criterion, (default, NaN: not applied)');
-			fielddisplay(self,'abstol','sea level change absolute convergence criterion, NaN: not applied');
+			fielddisplay(self,'reltol','sea level change relative convergence criterion (default, NaN: not applied)');
+			fielddisplay(self,'abstol','sea level change absolute convergence criterion(default, NaN: not applied)');
 			fielddisplay(self,'maxiter','maximum number of nonlinear iterations');
-			fielddisplay(self,'grdocean','does this planet have an ocean, if set to 1: global water mass is conserved in GRD module [default: 1]'); 
-			fielddisplay(self,'ocean_area_scaling','correction for model representation of ocean area [default: No correction]'); 
-			fielddisplay(self,'sealevelloading','enables surface loading from sea-level change (default 1)');
-			fielddisplay(self,'isgrd','compute GRD patterns (default 1)');
-			fielddisplay(self,'compute_bp_grd','compute GRD patterns for bottom pressure loads (default 1)');
+			fielddisplay(self,'grdocean','does this planet have an ocean, if set to 1: global water mass is conserved in GRD module (default: 1)'); 
+			fielddisplay(self,'ocean_area_scaling','correction for model representation of ocean area (default: No correction)'); 
+			fielddisplay(self,'sealevelloading','enables surface loading from sea-level change (default: 1)');
+			fielddisplay(self,'isgrd','compute GRD patterns (default: 1)');
+			fielddisplay(self,'compute_bp_grd','compute GRD patterns for bottom pressure loads (default: 1)');
 			fielddisplay(self,'runfrequency','how many time steps we skip before we run solidearthsettings solver during transient (default: 1)');
 			fielddisplay(self,'selfattraction','enables surface mass load to perturb the gravity field');
@@ -150,6 +151,6 @@
 			fielddisplay(self,'viscous','enables viscous deformation from surface loading');
 			fielddisplay(self,'rotation','enables polar motion to feedback on the GRD fields');
-			fielddisplay(self,'degacc','accuracy (default .01 deg) for numerical discretization of the Green''s functions');
-			fielddisplay(self,'timeacc','time accuracy (default 1 yr) for numerical discretization of the Green''s functions');
+			fielddisplay(self,'degacc','accuracy (default: .01 deg) for numerical discretization of the Green''s functions');
+			fielddisplay(self,'timeacc','time accuracy (default: 1 yr) for numerical discretization of the Green''s functions');
 			fielddisplay(self,'grdmodel','type of deformation model, 0 for no GRD, 1 for spherical GRD model (SESAW model), 2 for half-space planar GRD (visco-elastic model from Ivins)');
 			fielddisplay(self,'cross_section_shape','1: square-edged (default). 2: elliptical. See iedge in GiaDeflectionCore');
@@ -175,4 +176,6 @@
 			WriteData(fid,prefix,'object',self,'fieldname','cross_section_shape','name','md.solidearth.settings.cross_section_shape','format','Integer');
 		end % }}}
+		function self = extrude(self,md) % {{{
+		end % }}}
 		function savemodeljs(self,fid,modelname) % {{{
 		
@@ -191,6 +194,4 @@
 			writejsdouble(fid,[modelname '.solidearth.settings.cross_section_shape'],self.cross_section_shape);
 		end % }}}
-		function self = extrude(self,md) % {{{
-		end % }}}
 	end
 end
Index: /issm/trunk-jpl/src/m/classes/solidearthsettings.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/solidearthsettings.py	(revision 26298)
+++ /issm/trunk-jpl/src/m/classes/solidearthsettings.py	(revision 26299)
@@ -13,25 +13,25 @@
     """
 
-    def __init__(self, *args):  #{{{
+    def __init__(self, *args): #{{{
         self.reltol                 = 0
         self.abstol                 = 0
         self.maxiter                = 0
-        self.rigid                  = 0
-        self.elastic                = 0
-        self.rotation               = 0
+        self.selfattraction         = 1
+        self.elastic                = 1
+        self.viscous                = 1
+        self.rotation               = 1
+        self.grdocean               = 1
         self.ocean_area_scaling     = 0
         self.runfrequency           = 1 # How many time steps we skip before we run grd_core
-        self.computesealevelchange  = 1 # Will grd_core compute sea level?
-        self.isgrd                  = 1 # Will GRD patterns be computed?
-        self.compute_bp_grd         = 1 # Will GRD patterns for bottom pressures be computed?
+        self.sealevelloading        = 1 # Will sea-level loads be computed?
+        self.isgrd                  = 0 # Will GRD patterns be computed?
+        self.compute_bp_grd         = 0 # Will GRD patterns for bottom pressures be computed?
         self.degacc                 = 0 # Degree increment for resolution of Green tables
-        self.horiz                  = 0 # Compute horizontal displacement?
-        self.glfraction             = 1 # Barystatic contribution: full or fractional (default: fractional)
-        self.grdmodel               = 0 # GRD model (0 by default, 1 for elastic, 2 for Ivins)
+        self.timeacc                = 0 # Time step accuracy required to compute Green tables
+        self.horiz                  = 0 # Compute horizontal deformation?
+        self.grdmodel               = 0 # GRD model (0 by default, 1 for (visco-)elastic, 2 for Ivins)
         self.cross_section_shape    = 0 # Cross section only used when GRD model is Ivins
 
-        nargin = len(args)
-
-        if nargin == 0:
+        if len(args) == 0:
             self.setdefaultparameters()
         else:
@@ -39,24 +39,27 @@
     #}}}
 
-    def __repr__(self):  # {{{
+    def __repr__(self): #{{{
         s = '   solidearth settings:\n'
-        s += '{}\n'.format(fielddisplay(self, 'reltol', 'sea level rise relative convergence criterion, (default, NaN: not applied)'))
-        s += '{}\n'.format(fielddisplay(self, 'abstol', 'sea level rise absolute convergence criterion, NaN: not applied'))
+        s += '{}\n'.format(fielddisplay(self, 'reltol', 'sea level change relative convergence criterion (default, NaN: not applied)'))
+        s += '{}\n'.format(fielddisplay(self, 'abstol', 'sea level change absolute convergence criterion (default, NaN: not applied)'))
         s += '{}\n'.format(fielddisplay(self, 'maxiter', 'maximum number of nonlinear iterations'))
-        s += '{}\n'.format(fielddisplay(self, 'ocean_area_scaling', 'correction for model representation of ocean area [default: No correction]'))
-        s += '{}\n'.format(fielddisplay(self, 'computesealevelchange', 'compute sealevel change from GRD in addition to steric?) default 0'))
+        s += '{}\n'.format(fielddisplay(self, 'grdocean', 'does this planet have an ocean, if set to 1: global water mass is conserved in GRD module (default: 1)'))
+        s += '{}\n'.format(fielddisplay(self, 'ocean_area_scaling', 'correction for model representation of ocean area (default: No correction)'))
+        s += '{}\n'.format(fielddisplay(self, 'sealevelloading', 'enables surface loading from sea-level change (default: 1)'))
         s += '{}\n'.format(fielddisplay(self, 'isgrd', 'compute GRD patterns (default: 1'))
         s += '{}\n'.format(fielddisplay(self, 'compute_bp_grd', 'compute GRD patterns for bottom pressure loads (default 1)'))
         s += '{}\n'.format(fielddisplay(self, 'runfrequency', 'how many time steps we skip before we run solidearthsettings solver during transient (default: 1)'))
-        s += '{}\n'.format(fielddisplay(self, 'rigid', 'rigid earth graviational potential perturbation'))
-        s += '{}\n'.format(fielddisplay(self, 'elastic', 'elastic earth graviational potential perturbation'))
-        s += '{}\n'.format(fielddisplay(self, 'degacc', 'accuracy (default .01 deg) for numerical discretization of the Green\'s functions'))
-        s += '{}\n'.format(fielddisplay(self, 'glfraction', 'contribute fractionally (default, 1) to barystatic sea level'))
-        s += '{}\n'.format(fielddisplay(self, 'grdmodel', 'type of deformation model, 1 for elastic, 2 for visco-elastic from Ivins'))
+        s += '{}\n'.format(fielddisplay(self, 'selfattraction', 'enables surface mass load to perturb the gravity field'))
+        s += '{}\n'.format(fielddisplay(self, 'elastic', 'enables elastic deformation from surface loading'))
+        s += '{}\n'.format(fielddisplay(self, 'viscous', 'enables viscous deformation from surface loading'))
+        s += '{}\n'.format(fielddisplay(self, 'rotation', 'enables polar motion to feedback on the GRD fields'))
+        s += '{}\n'.format(fielddisplay(self, 'degacc', 'accuracy (default: .01 deg) for numerical discretization of the Green\'s functions'))
+        s += '{}\n'.format(fielddisplay(self, 'timeacc', 'time accuracy (default: 1 year) for numerical discretization of the Green\'s functions'))
+        s += '{}\n'.format(fielddisplay(self, 'grdmodel', 'type of deformation model, 0 for no GRD, 1 for spherical GRD model (SESAW model), 2 for half-space planar GRD (visco-elastic model from Ivins)'))
         s += '{}\n'.format(fielddisplay(self, 'cross_section_shape', '1: square-edged (default). 2: elliptical. See iedge in GiaDeflectionCore'))
         return s
     #}}}
 
-    def setdefaultparameters(self):  # {{{
+    def setdefaultparameters(self): #{{{
         # Convergence criterion: absolute, relative, and residual
         self.reltol = 0.01  # 1 percent
@@ -67,20 +70,20 @@
 
         # Computational flags
-        self.rigid = 1
+        self.selfattraction = 1
         self.elastic = 1
+        self.viscous = 1
         self.rotation = 1
+        self.grdocean = 1
         self.ocean_area_scaling = 0
-        self.compute_bp_grd = 1
+        self.compute_bp_grd = 0
         self.isgrd = 0
-        self.computesealevelchange = 1
+        self.sealevelloading = 1
 
         # Numerical discretization accuracy
         self.degacc = .01
+        self.timeacc = 1
 
         # How many time steps we skip before we run solidearthsettings solver during transient
         self.runfrequency = 1
-
-        # Fractional contribution
-        self.glfraction = 1
 
         # Horizontal displacement? (not on by default)
@@ -88,5 +91,5 @@
 
         # Cross section for Ivins model
-        self.cross_section_shape = 1 # Square as default (see iedde in GiaDeflectionCorex)
+        self.cross_section_shape = 1 # Square as default (see iedge in GiaDeflectionCorex)
 
         # No GRD model by default
@@ -94,7 +97,8 @@
     #}}}
 
-    def checkconsistency(self, md, solution, analyses):  # {{{
+    def checkconsistency(self, md, solution, analyses): #{{{
         if ('SealevelriseAnalysis' not in analyses) or (solution == 'TransientSolution' and not md.transient.isslc):
             return md
+
         md = checkfield(md, 'fieldname', 'solidearth.settings.reltol', 'size', [1])
         md = checkfield(md, 'fieldname', 'solidearth.settings.abstol', 'size', [1])
@@ -102,12 +106,13 @@
         md = checkfield(md, 'fieldname', 'solidearth.settings.runfrequency', 'size', [1], '>=', 1)
         md = checkfield(md, 'fieldname', 'solidearth.settings.degacc', 'size', [1], '>=', 1e-10)
+        md = checkfield(md, 'fieldname', 'solidearth.settings.timeacc', 'size', [1], '>', 0)
         md = checkfield(md, 'fieldname', 'solidearth.settings.horiz', 'NaN', 1, 'Inf', 1, 'values', [0, 1])
-        md = checkfield(md, 'fieldname', 'solidearth.settings.glfraction', 'values', [0, 1])
-        md = checkfield(md, 'fieldname', 'solidearth.settings.grdmodel', 'values', [1, 2])
+        md = checkfield(md, 'fieldname', 'solidearth.settings.grdmodel', '>=', 0, '<=', 2)
         md = checkfield(md, 'fieldname', 'solidearth.settings.cross_section_shape', 'numel', [1], 'values', [1, 2])
 
-        # Checks on computational flags
-        if self.elastic and not self.rigid:
-            raise Exception('solidearthsettings checkconsistency error message: need rigid on if elastic flag is set')
+        if self.elastic and not self.selfattraction:
+            raise Exception('solidearthsettings checkconsistency error message: need selfattraction on if elastic flag is set')
+        if self.viscous and not self.elastic:
+            raise Exception('solidearthsettings checkconsistency error message: need elastic on if viscous flag is set')
 
         # A GRD computation has been requested, make some checks on the nature of the meshes provided
@@ -119,4 +124,6 @@
                     if self.grdmodel == 1:
                         raise RuntimeException('model requires a 3D surface mesh to run GRD computations (change mesh from mesh2d to mesh3dsurface)')
+            if self.sealevelloading and not self.grdocean:
+                raise RuntimeException('solidearthsettings checkconsistency error message: need grdocean on if sealevelloading flag is set')
 
         if self.compute_bp_grd and not md.solidearth.settings.isgrd:
@@ -126,24 +133,26 @@
     #}}}
 
-    def marshall(self, prefix, md, fid):  #{{{
+    def marshall(self, prefix, md, fid): #{{{
         WriteData(fid, prefix, 'object', self, 'fieldname', 'reltol', 'name', 'md.solidearth.settings.reltol', 'format', 'Double')
         WriteData(fid, prefix, 'object', self, 'fieldname', 'abstol', 'name', 'md.solidearth.settings.abstol', 'format', 'Double')
         WriteData(fid, prefix, 'object', self, 'fieldname', 'maxiter', 'name', 'md.solidearth.settings.maxiter', 'format', 'Integer')
-        WriteData(fid, prefix, 'object', self, 'fieldname', 'rigid', 'name', 'md.solidearth.settings.rigid', 'format', 'Boolean')
+        WriteData(fid, prefix, 'object', self, 'fieldname', 'selfattraction', 'name', 'md.solidearth.settings.selfattraction', 'format', 'Boolean')
         WriteData(fid, prefix, 'object', self, 'fieldname', 'elastic', 'name', 'md.solidearth.settings.elastic', 'format', 'Boolean')
+        WriteData(fid, prefix, 'object', self, 'fieldname', 'viscous', 'name', 'md.solidearth.settings.viscous', 'format', 'Boolean')
         WriteData(fid, prefix, 'object', self, 'fieldname', 'rotation', 'name', 'md.solidearth.settings.rotation', 'format', 'Boolean')
-        WriteData(fid, prefix, 'object', self, 'fieldname', 'ocean_area_scaling', 'name', 'md.solidearth.settings.ocean_area_scaling', 'format', 'Boolean')
+        WriteData(fid, prefix, 'object', self, 'fieldname', 'grdocean', 'name', 'md.solidearth.settings.grdocean', 'format', 'Boolean')
+        WriteData(fid, prefix, 'object', self, 'fieldname', 'ocean_area_scaling', 'name', 'md.solidearth.settings.ocean_area_scaling', 'format', 'Integer')
         WriteData(fid, prefix, 'object', self, 'fieldname', 'runfrequency', 'name', 'md.solidearth.settings.runfrequency', 'format', 'Integer')
         WriteData(fid, prefix, 'object', self, 'fieldname', 'degacc', 'name', 'md.solidearth.settings.degacc', 'format', 'Double')
+        WriteData(fid, prefix, 'object', self, 'fieldname', 'timeacc', 'name', 'md.solidearth.settings.timeacc', 'format', 'Double', 'scale', md.constants.yts)
         WriteData(fid, prefix, 'object', self, 'fieldname', 'horiz', 'name', 'md.solidearth.settings.horiz', 'format', 'Integer')
-        WriteData(fid, prefix, 'object', self, 'fieldname', 'computesealevelchange', 'name', 'md.solidearth.settings.computesealevelchange', 'format', 'Integer')
+        WriteData(fid, prefix, 'object', self, 'fieldname', 'sealevelloading', 'name', 'md.solidearth.settings.sealevelloading', 'format', 'Integer')
         WriteData(fid, prefix, 'object', self, 'fieldname','isgrd', 'name', 'md.solidearth.settings.isgrd', 'format', 'Integer')
         WriteData(fid, prefix, 'object', self, 'fieldname', 'compute_bp_grd', 'name', 'md.solidearth.settings.compute_bp_grd', 'format', 'Integer')
-        WriteData(fid, prefix, 'object', self, 'fieldname','glfraction', 'name', 'md.solidearth.settings.glfraction', 'format', 'Integer')
         WriteData(fid, prefix, 'object', self, 'fieldname', 'grdmodel', 'name', 'md.solidearth.settings.grdmodel', 'format', 'Integer')
         WriteData(fid, prefix, 'object', self, 'fieldname', 'cross_section_shape', 'name', 'md.solidearth.settings.cross_section_shape', 'format', 'Integer')
     #}}}
 
-    def extrude(self, md):  #{{{
+    def extrude(self, md): #{{{
         return self
     #}}}
Index: /issm/trunk-jpl/test/NightlyRun/test2001.m
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test2001.m	(revision 26298)
+++ /issm/trunk-jpl/test/NightlyRun/test2001.m	(revision 26299)
@@ -1,4 +1,4 @@
 %Test Name: SquareSheetConstrainedGia2d
-%GIA test, inspired on test101. Running default GIA Ivins class.
+%GIA test, based off of test101. Running default GIA Ivins class.
 md=triangle(model(),'../Exp/Square.exp',100000.);
 md=setmask(md,'','');
Index: /issm/trunk-jpl/test/NightlyRun/test2001.py
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test2001.py	(revision 26298)
+++ /issm/trunk-jpl/test/NightlyRun/test2001.py	(revision 26299)
@@ -1,3 +1,4 @@
 #Test Name: SquareSheetConstrainedGia2d
+#GIA test, based off of test101. Running default GIA Ivins class.
 from socket import gethostname
 
@@ -19,4 +20,5 @@
 # GIA Ivins, 2 layer model
 md.solidearth.settings.grdmodel = 2
+md.solidearth.settings.isgrd = 1
 
 md.materials = materials('litho','ice')
@@ -28,6 +30,7 @@
 md.initialization.sealevel = np.zeros(md.mesh.numberofvertices)
 md.solidearth.settings.cross_section_shape = 1 # for square-edged x-section 
-md.solidearth.settings.computesealevelchange = 0 # do not compute sea level, only deformation
-md.solidearth.requested_outputs = ['Sealevel', 'SealevelUGrd']
+md.solidearth.settings.grdocean = 0 # do not compute sea level, only deformation
+md.solidearth.settings.sealevelloading = 0 # do not compute sea level, only deformation
+md.solidearth.requested_outputs = ['Sealevel', 'BedGRD']
 
 # Loading history
@@ -36,5 +39,5 @@
 # To get rid of default final_time, make sure final_time > start_time
 md.timestepping.final_time = 2400000  # 2,400 kyr
-md.masstransport.spcthickness np.array([
+md.masstransport.spcthickness = np.array([
     np.append(md.geometry.thickness, 0),
     np.append(md.geometry.thickness, 2400000)
@@ -50,9 +53,8 @@
 md.transient.isstressbalance = 0
 md.transient.isthermal = 0
-md.transient.ismasstransport = 0
-md.transient.isslc = 0
+md.transient.ismasstransport = 1
+md.transient.isslc = 1
 
 #Solve for GIA deflection
-md.cluster = generic('name', gethostname(), 'np', 1)
 md.cluster = generic('name', gethostname(), 'np', 3)
 md.verbose = verbose('11111111111')
@@ -63,3 +65,3 @@
 field_names = ['UGrd']
 field_tolerances = [1e-13]
-field_values = [md.results.TransientSolution[0].SealevelUGrd]
+field_values = [md.results.TransientSolution[1].BedGRD]