Index: /issm/trunk-jpl/src/m/classes/solidearth.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/solidearth.py	(revision 25763)
+++ /issm/trunk-jpl/src/m/classes/solidearth.py	(revision 25764)
@@ -20,6 +20,7 @@
 
     def __init__(self, *args):  #{{{
-        self.sealevel           = np.nan
+        self.initialsealevel    = np.nan
         self.settings           = solidearthsettings()
+        self.external           = []
         self.surfaceload        = surfaceload()
         self.lovenumbers        = lovenumbers()
@@ -41,12 +42,18 @@
     def __repr__(self):  # {{{
         s = '   solidearthinputs, forcings and settings:\n'
-        s += '{}\n'.format(fielddisplay(self, 'sealevel', 'current sea level (prior to computation) [m]'))
+        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'))
         s += '{}\n'.format(fielddisplay(self, 'requested_outputs', 'additional outputs requested'))
-        s += '{}\n'.format(self.settings)
-        s += '{}\n'.format(self.surfaceload)
-        s += '{}\n'.format(self.lovenumbers)
-        s += '{}\n'.format(self.rotational)
+        s += '{}\n'.format(fielddisplay(self, 'partitionice', 'ice partition vector for barystatic contribution'))
+        s += '{}\n'.format(fielddisplay(self, 'partitionhydro', 'hydro 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)
+        if self.external:
+            print(self.external)
         return s
     #}}}
@@ -63,4 +70,7 @@
         self.partitionhydro = []
 
+        # No external solutions by defalt
+        self.external = []
+
         # Earth radius
         self.planetradius = planetradius('earth')
@@ -70,5 +80,5 @@
         if ('SealevelriseAnalysis' not in analyses) or (solution == 'TransientSolution' and not md.transient.isslr):
             return md
-        md = checkfield(md, 'fieldname', 'solidearth.sealevel', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
+        md = checkfield(md, 'fieldname', 'solidearth.initialsealevel', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
         md = checkfield(md, 'fieldname', 'solidearth.requested_outputs', 'stringrow', 1)
 
@@ -77,4 +87,9 @@
         self.lovenumbers.checkconsistency(md, solution, analyses)
         self.rotational.checkconsistency(md, solution, analyses)
+        if self.external:
+            if not isinstance(self.external,'solidearthsolution'):
+                raise Exception('solidearth consistency check: external field should be a solidearthsolution')
+            end
+            self.external.checkconsistency(md,solution,analyses)
         return md
     #}}}
@@ -85,5 +100,5 @@
 
     def marshall(self, prefix, md, fid):  #{{{
-        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', 'initialsealevel', 'mattype', 1, 'format', 'DoubleMat', 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts)
         WriteData(fid, prefix, 'object', self, 'fieldname', 'planetradius', 'format', 'Double')
         WriteData(fid, prefix, 'object', self, 'fieldname', 'transitions', 'format', 'MatArray')
@@ -108,4 +123,6 @@
         self.lovenumbers.marshall(prefix, md, fid)
         self.rotational.marshall(prefix, md, fid)
+        if self.external:
+            self.external.marshall(prefix, md, fid)
 
         #process requested outputs
@@ -119,5 +136,5 @@
 
     def extrude(self, md): #{{{
-        self.sealevel = project3d(md, 'vector', self.sealevel, 'type', 'node')
+        self.initialsealevel = project3d(md, 'vector', self.initialsealevel, 'type', 'node')
         return self
     #}}}