Index: /issm/trunk-jpl/src/m/classes/giamme.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/giamme.py	(revision 25130)
+++ /issm/trunk-jpl/src/m/classes/giamme.py	(revision 25130)
@@ -0,0 +1,75 @@
+import numpy as np
+
+from checkfield import checkfield
+from fielddisplay import fielddisplay
+from project3d import project3d
+from WriteData import WriteData
+
+
+class giamme(object):
+    '''
+    GIAMME class definition
+
+        Usage:
+            gia = giamme() #gia class based on a multi-model ensemble (ex: Caron et al 2017 statistics)
+    '''
+
+    def __init__(self, *args): #{{{
+        self.modelid    = 0 #index into the multi-model ensemble
+        self.Ngia       = [] #geoid ensemble
+        self.Ugia       = [] #uplift ensemble
+        
+        nargin = len(args)
+
+        if nargin == 0:
+            self.setdefaultparameters()
+        else:
+            raise Exception('constructor not supported')
+    #}}}
+
+    def __repr__(self): # {{{
+        s = '   gia mme parameters:\n'
+        s += '{}\n'.format(fielddisplay(self, 'modelid', 'index into the multi-model ensemble, determines which field will be used.'))
+        s += '{}\n'.format(fielddisplay(self, 'Ngia', 'geoid (mm/yr).'))
+        s += '{}\n'.format(fielddisplay(self, 'Ugia', 'uplift (mm/yr).'))
+
+        return s
+    #}}}
+
+    def setdefaultparameters(self): # {{{
+        return self
+    #}}}
+
+    def checkconsistency(self, md, solution, analyses): # {{{
+        if not 'SealevelriseAnalysis' in analyses or solution == 'TransientSolution' and md.transient.istidallove == 0:
+            return md
+
+        md = checkfield(md, 'field', self.Ngia, 'NaN', 1, 'Inf', 1)
+        md = checkfield(md, 'field', self.Ugia, 'NaN', 1, 'Inf', 1)
+        md = checkfield(md, 'field', self.modelid, 'NaN', 1, 'Inf', 1, '>=', 0, '<=', len(self.Ngia))
+
+        return md
+    #}}}
+
+    def marshall(self, prefix, md, fid): #{{{
+        WriteData(fid, prefix, 'name', 'md.gia.model', 'data', 3, 'format', 'Integer')
+
+        if isinstance(self.Ngia, list) or self.Ngia.ndim == 1: #list or 1D numpy.ndarray
+            WriteData(fid, prefix, 'data', np.zeros(md.mesh.numberofvertices), 'format', 'DoubleMat', 'mattype', 1, 'name', 'md.gia.Ngia')
+            WriteData(fid, prefix, 'data', np.zeros(md.mesh.numberofvertices), 'format', 'DoubleMat', 'mattype', 1, 'name', 'md.gia.Ugia')
+            WriteData(fid, prefix, 'data', 1, 'name', 'md.gia.modelid', 'format', 'Double')
+            WriteData(fid, prefix, 'name', 'md.gia.nummodels', 'data', 1, 'format', 'Integer')
+        else: #2D numpy.ndarray
+            WriteData(fid, prefix, 'object ', self, 'fieldname', 'Ngia', 'format', 'DoubleMat', 'scale', 1 / 1000 / md.constants.yts)
+            WriteData(fid, prefix, 'object ', self, 'fieldname', 'Ugia', 'format', 'DoubleMat', 'scale', 1 / 1000 / md.constants.yts)
+            WriteData(fid, prefix, 'object ', self, 'fieldname', 'modelid', 'format', 'Double')
+            WriteData(fid, prefix, 'name', 'md.gia.nummodels', 'data', self.Ngia.shape[1],'format', 'Integer') #NOTE: <<numpy.ndarray>>.shape[1] gives length of second dimension
+    #}}}
+
+    def extrude(self, md): #{{{
+        for i in range(len(self.Ngia)):
+            self.Ngia[i] = project3d(md, 'vector', self.Ngia[i], 'type', 'node', 'layer', 1)
+            self.Ugia[i] = project3d(md, 'vector', self.Ugia[i], 'type', 'node', 'layer', 1)
+
+        return self
+    #}}}