source: issm/oecreview/Archive/21724-22754/ISSM-22266-22267.diff@ 22755

Index: ../trunk-jpl/test/Par/SquareThermal.py
===================================================================
--- ../trunk-jpl/test/Par/SquareThermal.py      (revision 22266)
+++ ../trunk-jpl/test/Par/SquareThermal.py      (revision 22267)
@@ -25,6 +25,9 @@
 
 print "      creating temperatures"
 md.initialization.temperature=(273.-20.)*numpy.ones((md.mesh.numberofvertices))
+md.initialization.pressure=numpy.zeros((md.mesh.numberofvertices,))
+md.initialization.waterfraction=numpy.zeros((md.mesh.numberofvertices,))
+md.initialization.watercolumn=numpy.zeros((md.mesh.numberofvertices,))
 
 print "      creating flow law parameter"
 md.materials.rheology_B=paterson(md.initialization.temperature)
@@ -32,7 +35,9 @@
 
 print "      creating surface mass balance"
 md.smb.mass_balance=numpy.ones((md.mesh.numberofvertices))/md.constants.yts    #1m/a
-md.basalforcings.melting_rate=0.*numpy.ones((md.mesh.numberofvertices))/md.constants.yts    #1m/a
+#md.basalforcings.melting_rate=0.*numpy.ones((md.mesh.numberofvertices))/md.constants.yts    #1m/a
+md.basalforcings.groundedice_melting_rate=0.*numpy.ones((md.mesh.numberofvertices))/md.constants.yts    #1m/a
+md.basalforcings.floatingice_melting_rate=0.*numpy.ones((md.mesh.numberofvertices))/md.constants.yts    #1m/a
 
 #Deal with boundary conditions:
 
@@ -40,6 +45,6 @@
 md=SetMarineIceSheetBC(md,'../Exp/SquareFront.exp')
 
 print "      boundary conditions for thermal model"
-md.thermal.spctemperature=md.initialization.temperature
+md.thermal.spctemperature[:]=md.initialization.temperature
 md.basalforcings.geothermalflux=numpy.zeros((md.mesh.numberofvertices)) 
 md.basalforcings.geothermalflux[numpy.nonzero(md.mask.groundedice_levelset>0.)[0]]=1.*10**-3    #1 mW/m^2
Index: ../trunk-jpl/test/Par/ISMIPE.par
===================================================================
--- ../trunk-jpl/test/Par/ISMIPE.par    (revision 22266)
+++ ../trunk-jpl/test/Par/ISMIPE.par    (revision 22267)
@@ -13,6 +13,7 @@
        md.geometry.base(i)=(1.-coeff)*data(point1,2)+coeff*data(point2,2);
        md.geometry.surface(i)=(1.-coeff)*data(point1,3)+coeff*data(point2,3);
 end
+
 md.geometry.thickness=md.geometry.surface-md.geometry.base;
 md.geometry.thickness(find(~md.geometry.thickness))=0.01;
 md.geometry.base=md.geometry.surface-md.geometry.thickness;
Index: ../trunk-jpl/test/Par/ISMIPE.py
===================================================================
--- ../trunk-jpl/test/Par/ISMIPE.py     (revision 22266)
+++ ../trunk-jpl/test/Par/ISMIPE.py     (revision 22267)
@@ -12,9 +12,10 @@
        y=md.mesh.y[i]
        point1=numpy.floor(y/100.)
        point2=numpy.minimum(point1+1,50)
-       coeff=(y-(point1-1.)*100.)/100.
+       coeff=(y-(point1)*100.)/100.
        md.geometry.base[i]=(1.-coeff)*data[point1,1]+coeff*data[point2,1]
        md.geometry.surface[i]=(1.-coeff)*data[point1,2]+coeff*data[point2,2]
+
 md.geometry.thickness=md.geometry.surface-md.geometry.base
 md.geometry.thickness[numpy.nonzero(numpy.logical_not(md.geometry.thickness))]=0.01
 md.geometry.base=md.geometry.surface-md.geometry.thickness
Index: ../trunk-jpl/test/NightlyRun/test2425.py
===================================================================
--- ../trunk-jpl/test/NightlyRun/test2425.py    (nonexistent)
+++ ../trunk-jpl/test/NightlyRun/test2425.py    (revision 22267)
@@ -0,0 +1,50 @@
+#Test Name: SquareSheetShelfGroundingLine2dSoft
+import numpy as np
+from model import *
+from socket import gethostname
+from triangle import *
+from setmask import *
+from parameterize import *
+from setflowequation import *
+from solve import *
+
+md = triangle(model(),'../Exp/Square.exp',150000.)
+md = setmask(md,'../Exp/SquareShelf.exp','')
+md = parameterize(md,'../Par/SquareSheetShelf.py')
+md = setflowequation(md,'SSA','all')
+md.initialization.vx[:] = 0.
+md.initialization.vy[:] = 0.
+md.geometry.base = -700. - (md.mesh.y - 500000.) / 1000.
+md.geometry.bed = -700. - (md.mesh.y - 500000.) / 1000.
+md.geometry.thickness[:] = 1300.
+md.geometry.surface = md.geometry.base + md.geometry.thickness
+md.transient.isstressbalance = 1
+md.transient.isgroundingline = 1
+md.groundingline.migration = 'AggressiveMigration'
+
+md.timestepping.time_step = .1
+md.timestepping.final_time = 1
+
+md.cluster = generic('name',gethostname(),'np',3)
+md = solve(md,'Transient')
+vel1 = md.results.TransientSolution[-1].Vel
+
+#get same results with offset in bed and sea level: 
+md.geometry.base = -700. - (md.mesh.y - 500000.) / 1000.
+md.geometry.bed  = -700. - (md.mesh.y - 500000.) / 1000.
+md.geometry.thickness[:] = 1300.
+md.geometry.surface = md.geometry.base + md.geometry.thickness
+
+md.geometry.base = md.geometry.base + 1000
+md.geometry.bed = md.geometry.bed + 1000
+md.geometry.surface = md.geometry.surface + 1000
+md.slr.sealevel = 1000 * np.ones((md.mesh.numberofvertices,))
+
+md = solve(md,'Transient','checkconsistency','no')
+vel2 = md.results.TransientSolution[-1].Vel
+
+#Fields and tolerances to track changes
+field_names = ['Vel','Veloffset']
+field_tolerances = [1e-13,1e-13]
+field_values = [vel1,vel2]
+
Index: ../trunk-jpl/test/NightlyRun/test342.py
===================================================================
--- ../trunk-jpl/test/NightlyRun/test342.py     (nonexistent)
+++ ../trunk-jpl/test/NightlyRun/test342.py     (revision 22267)
@@ -0,0 +1,35 @@
+#Test Name: SquareSheetTherSteaPlume
+import numpy as np
+from model import *
+from socket import gethostname
+from triangle import *
+from setmask import *
+from parameterize import *
+from setflowequation import *
+from solve import *
+from plumebasalforcings import *
+
+md = triangle(model(),'../Exp/Square.exp',180000.)
+md = setmask(md,'','')
+md = parameterize(md,'../Par/SquareSheetConstrained.py')
+md.basalforcings = plumebasalforcings()
+md.basalforcings = md.basalforcings.setdefaultparameters()
+md.basalforcings.floatingice_melting_rate = np.zeros((md.mesh.numberofvertices,))
+md.basalforcings.groundedice_melting_rate = np.zeros((md.mesh.numberofvertices,))
+md.basalforcings.plumex = 500000
+md.basalforcings.plumey = 500000
+md.extrude(3,1.)
+md = setflowequation(md,'SSA','all')
+md.timestepping.time_step = 0.
+md.thermal.requested_outputs = ['default','BasalforcingsGeothermalflux']
+md.cluster = generic('name',gethostname(),'np',3)
+md = solve(md,'Thermal')
+
+#Fields and tolerances to track changes
+field_names      = ['Temperature','BasalforcingsGroundediceMeltingRate','BasalforcingsGeothermalflux']
+field_tolerances = [1e-13,1e-8,1e-13]
+field_values = [
+       md.results.ThermalSolution.Temperature,
+       md.results.ThermalSolution.BasalforcingsGroundediceMeltingRate,
+       md.results.ThermalSolution.BasalforcingsGeothermalflux,
+       ]
Index: ../trunk-jpl/test/NightlyRun/test261.py
===================================================================
--- ../trunk-jpl/test/NightlyRun/test261.py     (nonexistent)
+++ ../trunk-jpl/test/NightlyRun/test261.py     (revision 22267)
@@ -0,0 +1,38 @@
+#Test Name: SquareShelfConstrainedTranEnhanced
+import numpy as np
+from model import *
+from socket import gethostname
+from triangle import *
+from setmask import *
+from parameterize import *
+from setflowequation import *
+from solve import *
+from matenhancedice import *
+
+md = triangle(model(),'../Exp/Square.exp',180000.)
+md = setmask(md,'all','')
+md = parameterize(md,'../Par/SquareShelfConstrained.py')
+md = md.extrude(3,1.)
+md = setflowequation(md,'SSA','all')
+md.cluster = generic('name',gethostname(),'np',3)
+md.materials = matenhancedice()
+md.materials.rheology_B = 3.15e8 * np.ones(md.mesh.numberofvertices,)
+md.materials.rheology_n = 3 * np.ones(md.mesh.numberofelements,)
+md.materials.rheology_E = np.ones(md.mesh.numberofvertices,)
+md.transient.isstressbalance = 1
+md.transient.ismasstransport = 0
+md.transient.issmb = 1
+md.transient.isthermal = 1
+md.transient.isgroundingline = 0
+md = solve(md,'Transient')
+
+#Fields and tolerances to track changes
+field_names      = ['Vx','Vy','Vel','Temperature','BasalforcingsGroundediceMeltingRate']
+field_tolerances = [1e-13,1e-13,1e-13,1e-13,1e-13]
+field_values = [
+       md.results.TransientSolution[0].Vx,
+       md.results.TransientSolution[0].Vy,
+       md.results.TransientSolution[0].Vel,
+       md.results.TransientSolution[0].Temperature,
+       md.results.TransientSolution[0].BasalforcingsGroundediceMeltingRate,
+       ]
Index: ../trunk-jpl/test/NightlyRun/test441.py
===================================================================
--- ../trunk-jpl/test/NightlyRun/test441.py     (nonexistent)
+++ ../trunk-jpl/test/NightlyRun/test441.py     (revision 22267)
@@ -0,0 +1,93 @@
+#Test Name: MISMIP3D
+import numpy as np
+from model import *
+from socket import gethostname
+from triangle import *
+from setmask import *
+from parameterize import *
+from setflowequation import *
+from solve import *
+
+md = triangle(model(),'../Exp/Square.exp',100000.)
+md = setmask(md,'../Exp/SquareShelf.exp','')
+md = parameterize(md,'../Par/SquareSheetShelf.py')
+md.initialization.vx[:] = 1.
+md.initialization.vy[:] = 1.
+md.geometry.thickness[:] = 500. - md.mesh.x / 10000.
+md.geometry.bed = -100. - md.mesh.x / 1000.
+md.geometry.base = -md.geometry.thickness * md.materials.rho_ice / md.materials.rho_water
+md.mask.groundedice_levelset = md.geometry.thickness + md.materials.rho_water / md.materials.rho_ice * md.geometry.bed
+pos = np.array(np.where(md.mask.groundedice_levelset >= 0.))
+md.geometry.base[pos] = md.geometry.bed[pos]
+md.geometry.surface = md.geometry.base + md.geometry.thickness
+md = setflowequation(md,'SSA','all')
+
+#Boundary conditions:
+md.mask.ice_levelset = -np.ones((md.mesh.numberofvertices,))
+md.mask.ice_levelset[np.where(md.mesh.x == max(md.mesh.x))] = 0.
+md.stressbalance.spcvx[:] = float('Nan')
+md.stressbalance.spcvy[:] = float('Nan')
+md.stressbalance.spcvz[:] = float('Nan')
+posA = np.intersect1d(np.array(np.where(md.mesh.y < 1000000.1)),np.array(np.where(md.mesh.y > 999999.9)))
+posB = np.intersect1d(np.array(np.where(md.mesh.y < 0.1)),np.array(np.where(md.mesh.y > -0.1)))
+pos = np.unique(np.concatenate((posA,posB)))
+md.stressbalance.spcvy[pos] = 0.
+pos2 = np.intersect1d(np.array(np.where(md.mesh.x < 0.1)), np.array(np.where(md.mesh.x > -0.1)))
+md.stressbalance.spcvx[pos2] = 0.
+md.stressbalance.spcvy[pos2] = 0.
+
+md.materials.rheology_B = 1. / ((10**-25)**(1./3.)) * np.ones((md.mesh.numberofvertices,))
+md.materials.rheology_law = 'None'
+md.friction.coefficient[:] = np.sqrt(1e7) * np.ones((md.mesh.numberofvertices,))
+md.friction.p = 3. * np.ones((md.mesh.numberofelements,))
+md.smb.mass_balance[:] = 1.
+md.basalforcings.groundedice_melting_rate[:] = 0.
+md.basalforcings.floatingice_melting_rate[:] = 30.
+md.transient.isthermal = 0
+md.transient.isstressbalance = 1
+md.transient.isgroundingline = 1
+md.transient.ismasstransport = 1
+md.transient.issmb = 1
+md.transient.requested_outputs = ['default','BasalforcingsFloatingiceMeltingRate']
+md.groundingline.migration = 'SubelementMigration2'
+md.timestepping.final_time = 30.
+md.timestepping.time_step = 10.
+
+md.cluster = generic('name',gethostname(),'np',3)
+md = solve(md,'Transient')
+
+#Fields and tolerances to track changes
+field_names     = [
+       'Bed1','Surface1','Thickness1','Floatingice1','Vx1','Vy1','Pressure1','FloatingiceMeltingrate1',
+       'Bed2','Surface2','Thickness2','Floatingice2','Vx2','Vy2','Pressure2','FloatingiceMeltingrate2',
+       'Bed3','Surface3','Thickness3','Floatingice3','Vx3','Vy3','Pressure3','FloatingiceMeltingrate3']
+field_tolerances = [
+       2e-11,5e-12,2e-11,1e-11,5e-10,1e-08,1e-13,1e-13,
+       3e-11,3e-11,9e-10,7e-11,1e-09,5e-08,1e-10,1e-13,
+       1e-10,3e-11,1e-10,7e-11,1e-09,5e-08,1e-10,1e-13]
+field_values = [
+       md.results.TransientSolution[0].Base,
+       md.results.TransientSolution[0].Surface,
+       md.results.TransientSolution[0].Thickness,
+       md.results.TransientSolution[0].MaskGroundediceLevelset,
+       md.results.TransientSolution[0].Vx,
+       md.results.TransientSolution[0].Vy,
+       md.results.TransientSolution[0].Pressure,
+       md.results.TransientSolution[0].BasalforcingsFloatingiceMeltingRate,
+       md.results.TransientSolution[1].Base,
+       md.results.TransientSolution[1].Surface,
+       md.results.TransientSolution[1].Thickness,
+       md.results.TransientSolution[1].MaskGroundediceLevelset,
+       md.results.TransientSolution[1].Vx,
+       md.results.TransientSolution[1].Vy,
+       md.results.TransientSolution[1].Pressure,
+       md.results.TransientSolution[1].BasalforcingsFloatingiceMeltingRate,
+       md.results.TransientSolution[2].Base,
+       md.results.TransientSolution[2].Surface,
+       md.results.TransientSolution[2].Thickness,
+       md.results.TransientSolution[2].MaskGroundediceLevelset,
+       md.results.TransientSolution[2].Vx,
+       md.results.TransientSolution[2].Vy,
+       md.results.TransientSolution[2].Pressure,
+       md.results.TransientSolution[2].BasalforcingsFloatingiceMeltingRate,
+       ]
Index: ../trunk-jpl/test/NightlyRun/test343.py
===================================================================
--- ../trunk-jpl/test/NightlyRun/test343.py     (nonexistent)
+++ ../trunk-jpl/test/NightlyRun/test343.py     (revision 22267)
@@ -0,0 +1,66 @@
+#Test Name: SquareSheetConstrainedSmbGradientsEla2d
+import numpy as np
+from model import *
+from socket import gethostname
+from triangle import *
+from setmask import *
+from parameterize import *
+from setflowequation import *
+from solve import *
+from SMBgradientsela import *
+
+md = triangle(model(),'../Exp/Square.exp',150000.)
+md = setmask(md,'','')
+md = parameterize(md,'../Par/SquareSheetConstrained.py')
+md = setflowequation(md,'SSA','all')
+md.smb = SMBgradientsela()
+md.smb.ela = 1500. * np.ones((md.mesh.numberofvertices+1,))
+md.smb.b_pos = 0.002 * np.ones((md.mesh.numberofvertices+1,))
+md.smb.b_neg = 0.005 * np.ones((md.mesh.numberofvertices+1,))
+md.smb.b_max = 4. * (md.materials.rho_freshwater / md.materials.rho_ice) * np.ones((md.mesh.numberofvertices+1,))
+md.smb.b_min = -4. * (md.materials.rho_freshwater / md.materials.rho_ice) * np.ones((md.mesh.numberofvertices+1,))
+
+#Change geometry
+md.geometry.thickness = md.geometry.surface * 30.
+md.geometry.surface = md.geometry.base + md.geometry.thickness
+
+#Transient options
+md.transient.requested_outputs = ['default','TotalSmb']
+md.cluster = generic('name',gethostname(),'np',3)
+md = solve(md,'Transient')
+
+#Fields and tolerances to track changes
+field_names     = [
+       'Vx1','Vy1','Vel1','Bed1','Surface1','Thickness1','SMB1','TotalSmb1',
+       'Vx2','Vy2','Vel2','Bed2','Surface2','Thickness2','SMB2','TotalSmb2',
+       'Vx3','Vy3','Vel3','Bed3','Surface3','Thickness3','SMB3','TotalSmb3']
+field_tolerances = [
+       1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,
+       1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,
+       1e-12,1e-12,1e-12,1e-13,1e-13,1e-13,1.5e-13,1e-13]
+field_values = [
+       md.results.TransientSolution[0].Vx,
+       md.results.TransientSolution[0].Vy,
+       md.results.TransientSolution[0].Vel,
+       md.results.TransientSolution[0].Base,
+       md.results.TransientSolution[0].Surface,
+       md.results.TransientSolution[0].Thickness,
+       md.results.TransientSolution[0].SmbMassBalance,
+       md.results.TransientSolution[0].TotalSmb,
+       md.results.TransientSolution[1].Vx,
+       md.results.TransientSolution[1].Vy,
+       md.results.TransientSolution[1].Vel,
+       md.results.TransientSolution[1].Base,
+       md.results.TransientSolution[1].Surface,
+       md.results.TransientSolution[1].Thickness,
+       md.results.TransientSolution[1].TotalSmb,
+       md.results.TransientSolution[1].SmbMassBalance,
+       md.results.TransientSolution[2].Vx,
+       md.results.TransientSolution[2].Vy,
+       md.results.TransientSolution[2].Vel,
+       md.results.TransientSolution[2].Base,
+       md.results.TransientSolution[2].Surface,
+       md.results.TransientSolution[2].Thickness,
+       md.results.TransientSolution[2].SmbMassBalance,
+       md.results.TransientSolution[2].TotalSmb
+       ]
Index: ../trunk-jpl/test/NightlyRun/test540.py
===================================================================
--- ../trunk-jpl/test/NightlyRun/test540.py     (nonexistent)
+++ ../trunk-jpl/test/NightlyRun/test540.py     (revision 22267)
@@ -0,0 +1,66 @@
+#Test Name: PigTranCalvingDevSSA2d
+from model import *
+from socket import gethostname
+from triangle import *
+from setmask import *
+from parameterize import *
+from setflowequation import *
+from solve import *
+from calvingvonmises import *
+
+md = triangle(model(),'../Exp/Pig.exp',10000.)
+md = setmask(md,'../Exp/PigShelves.exp','../Exp/PigIslands.exp')
+md = parameterize(md,'../Par/Pig.py')
+md = setflowequation(md,'SSA','all')
+md.timestepping.time_step = 2
+md.timestepping.final_time = 50
+
+#calving parameters
+md.mask.ice_levelset = 1e4 * (md.mask.ice_levelset + 0.5)
+md.calving = calvingvonmises()
+md.calving.meltingrate = np.zeros((md.mesh.numberofvertices,))
+md.transient.ismovingfront = 1
+md.levelset.spclevelset = float('NaN') * np.ones((md.mesh.numberofvertices,))
+pos = np.where(md.mesh.vertexonboundary)
+md.levelset.spclevelset[pos] = md.mask.ice_levelset[pos]
+
+md.cluster = generic('name',gethostname(),'np',2)
+md = solve(md,'Transient')
+
+#Fields and tolerances to track changes
+field_names = [
+       'Vx1' ,'Vy1' ,'Vel1' ,'Pressure1' ,'Bed1' ,'Surface1' ,'Thickness1' ,'MaskIceLevelset1' ,
+       'Vx2' ,'Vy2' ,'Vel2' ,'Pressure2' ,'Bed2' ,'Surface2' ,'Thickness2' ,'MaskIceLevelset2' ,
+       'Vx10','Vy10','Vel10','Pressure10','Bed10','Surface10','Thickness10','MaskIceLevelset10',
+       ]
+field_tolerances = [
+       1e-12,2e-12,2e-12,1e-13,1e-13,1e-13,1e-13,1e-13,
+       1e-12,1e-12,1e-12,1e-13,1e-13,1e-13,1e-13,1e-12,
+       1e-11,1e-11,1e-11,1e-11,1e-11,1e-11,1e-11,1e-9,
+       ]
+field_values = [
+       md.results.TransientSolution[0].Vx,
+       md.results.TransientSolution[0].Vy,
+       md.results.TransientSolution[0].Vel,
+       md.results.TransientSolution[0].Pressure,
+       md.results.TransientSolution[0].Base,
+       md.results.TransientSolution[0].Surface,
+       md.results.TransientSolution[0].Thickness,
+       md.results.TransientSolution[0].MaskIceLevelset,
+       md.results.TransientSolution[1].Vx,
+       md.results.TransientSolution[1].Vy,
+       md.results.TransientSolution[1].Vel,
+       md.results.TransientSolution[1].Pressure,
+       md.results.TransientSolution[1].Base,
+       md.results.TransientSolution[1].Surface,
+       md.results.TransientSolution[1].Thickness,
+       md.results.TransientSolution[1].MaskIceLevelset,
+       md.results.TransientSolution[9].Vx,
+       md.results.TransientSolution[9].Vy,
+       md.results.TransientSolution[9].Vel,
+       md.results.TransientSolution[9].Pressure,
+       md.results.TransientSolution[9].Base,
+       md.results.TransientSolution[9].Surface,
+       md.results.TransientSolution[9].Thickness,
+       md.results.TransientSolution[9].MaskIceLevelset,
+       ]
Index: ../trunk-jpl/test/NightlyRun/test442.py
===================================================================
--- ../trunk-jpl/test/NightlyRun/test442.py     (nonexistent)
+++ ../trunk-jpl/test/NightlyRun/test442.py     (revision 22267)
@@ -0,0 +1,97 @@
+#Test Name: MISMIP3DHO
+import numpy as np
+from model import *
+from socket import gethostname
+from triangle import *
+from setmask import *
+from parameterize import *
+from setflowequation import *
+from solve import *
+
+md = triangle(model(),'../Exp/Square.exp',100000.)
+md = setmask(md,'../Exp/SquareShelf.exp','')
+md = parameterize(md,'../Par/SquareSheetShelf.py')
+md.initialization.vx[:] = 1.
+md.initialization.vy[:] = 1.
+md.geometry.thickness[:] = 500. - md.mesh.x / 10000.
+md.geometry.bed = -100. - md.mesh.x / 1000.
+md.geometry.base = -md.geometry.thickness * md.materials.rho_ice / md.materials.rho_water
+md.mask.groundedice_levelset = md.geometry.thickness + md.materials.rho_water / md.materials.rho_ice * md.geometry.bed
+pos = np.where(md.mask.groundedice_levelset >= 0.)
+md.geometry.base[pos] = md.geometry.bed[pos]
+md.geometry.surface = md.geometry.base + md.geometry.thickness
+md = md.extrude(4,1.)
+md = setflowequation(md,'HO','all')
+
+#Boundary conditions:
+md.mask.ice_levelset = -np.ones((md.mesh.numberofvertices,))
+md.mask.ice_levelset[np.where(md.mesh.x == max(md.mesh.x))] = 0.
+md.stressbalance.spcvx[:] = float('Nan')
+md.stressbalance.spcvy[:] = float('Nan')
+md.stressbalance.spcvz[:] = float('Nan')
+posA = np.intersect1d(np.array(np.where(md.mesh.y < 1000000.1)),np.array(np.where(md.mesh.y > 999999.9)))
+posB = np.intersect1d(np.array(np.where(md.mesh.y < 0.1)),np.array(np.where(md.mesh.y > -0.1)))
+pos = np.unique(np.concatenate((posA,posB)))
+md.stressbalance.spcvy[pos] = 0.
+pos2 = np.intersect1d(np.array(np.where(md.mesh.x < 0.1)), np.array(np.where(md.mesh.x > -0.1)))
+md.stressbalance.spcvx[pos2] = 0.
+md.stressbalance.spcvy[pos2] = 0.
+
+md.materials.rheology_B = 1. / ((10**-25)**(1./3.)) * np.ones((md.mesh.numberofvertices,))
+md.materials.rheology_law = 'None'
+md.friction.coefficient[:] = np.sqrt(1e7) * np.ones((md.mesh.numberofvertices,))
+md.friction.p = 3. * np.ones((md.mesh.numberofelements,))
+md.smb.mass_balance[:] = 1.
+md.basalforcings.groundedice_melting_rate[:] = 0.
+md.basalforcings.floatingice_melting_rate[:] = 30.
+md.transient.isthermal = 0
+md.transient.isstressbalance = 1
+md.transient.isgroundingline = 1
+md.transient.ismasstransport = 1
+md.transient.issmb = 1
+md.transient.requested_outputs = ['default','BasalforcingsFloatingiceMeltingRate']
+md.groundingline.migration = 'SubelementMigration2'
+md.timestepping.final_time = 30
+md.timestepping.time_step = 10
+
+md.cluster = generic('name',gethostname(),'np',3)
+md = solve(md,'Transient')
+
+#Fields and tolerances to track changes
+field_names     = [
+       'Bed1','Surface1','Thickness1','Floatingice1','Vx1','Vy1','Vz1','Pressure1','FloatingiceMeltingrate1',
+       'Bed2','Surface2','Thickness2','Floatingice2','Vx2','Vy2','Vz2','Pressure2','FloatingiceMeltingrate2',
+       'Bed3','Surface3','Thickness3','Floatingice3','Vx3','Vy3','Vz3','Pressure3','FloatingiceMeltingrate3']
+field_tolerances = [
+       2e-11,5e-12,2e-11,1e-11,5e-10,3e-08,6e-10,1e-13,1e-13,
+       3e-11,3e-11,9e-10,7e-11,6e-09,1e-07,1e-09,1e-10,1e-13,
+       1e-8,2e-08,7e-09,2e-7 ,1e-03,8e-04,2e-09,1e-10,1e-13]
+field_values = [
+       md.results.TransientSolution[0].Base,
+       md.results.TransientSolution[0].Surface,
+       md.results.TransientSolution[0].Thickness,
+       md.results.TransientSolution[0].MaskGroundediceLevelset,
+       md.results.TransientSolution[0].Vx,
+       md.results.TransientSolution[0].Vy,
+       md.results.TransientSolution[0].Vz,
+       md.results.TransientSolution[0].Pressure,
+       md.results.TransientSolution[0].BasalforcingsFloatingiceMeltingRate,
+       md.results.TransientSolution[1].Base,
+       md.results.TransientSolution[1].Surface,
+       md.results.TransientSolution[1].Thickness,
+       md.results.TransientSolution[1].MaskGroundediceLevelset,
+       md.results.TransientSolution[1].Vx,
+       md.results.TransientSolution[1].Vy,
+       md.results.TransientSolution[1].Vz,
+       md.results.TransientSolution[1].Pressure,
+       md.results.TransientSolution[1].BasalforcingsFloatingiceMeltingRate,
+       md.results.TransientSolution[2].Base,
+       md.results.TransientSolution[2].Surface,
+       md.results.TransientSolution[2].Thickness,
+       md.results.TransientSolution[2].MaskGroundediceLevelset,
+       md.results.TransientSolution[2].Vx,
+       md.results.TransientSolution[2].Vy,
+       md.results.TransientSolution[2].Vz,
+       md.results.TransientSolution[2].Pressure,
+       md.results.TransientSolution[2].BasalforcingsFloatingiceMeltingRate,
+       ]
Index: ../trunk-jpl/test/NightlyRun/test344.py
===================================================================
--- ../trunk-jpl/test/NightlyRun/test344.py     (nonexistent)
+++ ../trunk-jpl/test/NightlyRun/test344.py     (revision 22267)
@@ -0,0 +1,73 @@
+#Test Name: SquareSheetConstrainedSmbGradientsEla3d
+import numpy as np
+from model import *
+from socket import gethostname
+from triangle import *
+from setmask import *
+from parameterize import *
+from setflowequation import *
+from solve import *
+from SMBgradientsela import *
+
+md = triangle(model(),'../Exp/Square.exp',150000.)
+md = setmask(md,'','')
+md = parameterize(md,'../Par/SquareSheetConstrained.py')
+
+#Change geometry
+md.geometry.thickness = md.geometry.surface * 30.
+md.geometry.surface = md.geometry.base + md.geometry.thickness
+
+md = md.extrude(3,1.)
+md = setflowequation(md,'HO','all')
+md.smb = SMBgradientsela()
+md.smb.ela = 1500. * np.ones((md.mesh.numberofvertices+1,))
+md.smb.b_pos = 0.002 * np.ones((md.mesh.numberofvertices+1,))
+md.smb.b_neg = 0.005 * np.ones((md.mesh.numberofvertices+1,))
+md.smb.b_max = 4. * (md.materials.rho_freshwater / md.materials.rho_ice) * np.ones((md.mesh.numberofvertices+1,))
+md.smb.b_min = -4. * (md.materials.rho_freshwater / md.materials.rho_ice) * np.ones((md.mesh.numberofvertices+1,))
+
+
+#Transient options
+md.transient.requested_outputs = ['default','TotalSmb']
+md.cluster = generic('name',gethostname(),'np',3)
+md = solve(md,'Transient')
+
+#Fields and tolerances to track changes
+field_names     = ['Vx1','Vy1','Vz1','Vel1','Bed1','Surface1','Thickness1','Temperature1','SMB1','TotalSmb1',
+ 'Vx2','Vy2','Vz2','Vel2','Bed2','Surface2','Thickness2','Temperature2','SMB2','TotalSmb2',
+ 'Vx3','Vy3','Vz3','Vel3','Bed3','Surface3','Thickness3','Temperature3','SMB3','TotalSmb3']
+field_tolerances = [1e-09,1e-09,1e-09,1e-09,1e-10,1e-10,1e-10,1e-10,1e-10,1e-10,
+       1e-09,1e-09,1e-10,1e-09,1e-10,1e-10,1e-10,1e-10,1e-10,1e-10,
+       1e-09,5e-09,1e-09,1e-09,1e-10,1e-10,1e-10,1e-10,1e-10,1e-10]
+field_values = [
+       md.results.TransientSolution[0].Vx,
+       md.results.TransientSolution[0].Vy,
+       md.results.TransientSolution[0].Vz,
+       md.results.TransientSolution[0].Vel,
+       md.results.TransientSolution[0].Base,
+       md.results.TransientSolution[0].Surface,
+       md.results.TransientSolution[0].Thickness,
+       md.results.TransientSolution[0].Temperature,
+       md.results.TransientSolution[0].SmbMassBalance,
+       md.results.TransientSolution[0].TotalSmb,
+       md.results.TransientSolution[1].Vx,
+       md.results.TransientSolution[1].Vy,
+       md.results.TransientSolution[1].Vz,
+       md.results.TransientSolution[1].Vel,
+       md.results.TransientSolution[1].Base,
+       md.results.TransientSolution[1].Surface,
+       md.results.TransientSolution[1].Thickness,
+       md.results.TransientSolution[1].Temperature,
+       md.results.TransientSolution[1].SmbMassBalance,
+       md.results.TransientSolution[1].TotalSmb,
+       md.results.TransientSolution[2].Vx,
+       md.results.TransientSolution[2].Vy,
+       md.results.TransientSolution[2].Vz,
+       md.results.TransientSolution[2].Vel,
+       md.results.TransientSolution[2].Base,
+       md.results.TransientSolution[2].Surface,
+       md.results.TransientSolution[2].Thickness,
+       md.results.TransientSolution[2].Temperature,
+       md.results.TransientSolution[2].SmbMassBalance,
+       md.results.TransientSolution[2].TotalSmb,
+       ]
Index: ../trunk-jpl/test/NightlyRun/test460.py
===================================================================
--- ../trunk-jpl/test/NightlyRun/test460.py     (nonexistent)
+++ ../trunk-jpl/test/NightlyRun/test460.py     (revision 22267)
@@ -0,0 +1,38 @@
+#Test Name: SquareSheetShelfStressFSEstar
+import numpy as np
+from model import *
+from socket import gethostname
+from triangle import *
+from setmask import *
+from parameterize import *
+from setflowequation import *
+from solve import *
+from matestar import *
+
+md=triangle(model(),'../Exp/Square.exp',180000.)
+md=setmask(md,'../Exp/SquareShelf.exp','')
+md=parameterize(md,'../Par/SquareSheetShelf.py')
+md = md.extrude(3,1.)
+md.materials = matestar()
+md.materials.rheology_B = 3.15e8 * np.ones((md.mesh.numberofvertices,))
+md.materials.rheology_Ec = np.ones((md.mesh.numberofvertices,))
+md.materials.rheology_Es = 3 * np.ones((md.mesh.numberofvertices,))
+md.cluster = generic('name',gethostname(),'np',3)
+
+#Go solve
+field_names=[]
+field_tolerances=[]
+field_values=[]
+#md.initialization.pressure=md.constants.g*md.materials.rho_ice*(md.geometry.surface-md.mesh.y);
+for i in ['SSA','HO','FS']:
+       md = setflowequation(md,i,'all')
+       md = solve(md,'Stressbalance')
+       field_names     = field_names + ['Vx'+i,'Vy'+i,'Vz'+i,'Vel'+i,'Pressure'+i]
+       field_tolerances = field_tolerances + [6e-07,6e-07,2e-06,1e-06,5e-07]
+       field_values = field_values + [
+                       md.results.StressbalanceSolution.Vx,
+                       md.results.StressbalanceSolution.Vy,
+                       md.results.StressbalanceSolution.Vz,
+                       md.results.StressbalanceSolution.Vel,
+                       md.results.StressbalanceSolution.Pressure,
+                       ]
Index: ../trunk-jpl/test/NightlyRun/test461.py
===================================================================
--- ../trunk-jpl/test/NightlyRun/test461.py     (nonexistent)
+++ ../trunk-jpl/test/NightlyRun/test461.py     (revision 22267)
@@ -0,0 +1,53 @@
+#Test Name: SquareSheetShelfThermalFSEstar
+import numpy as np
+from model import *
+from socket import gethostname
+from triangle import *
+from setmask import *
+from parameterize import *
+from setflowequation import *
+from solve import *
+from matestar import *
+
+md = triangle(model(),'../Exp/Square.exp',180000.)
+md = setmask(md,'../Exp/SquareShelf.exp','')
+md = parameterize(md,'../Par/SquareSheetShelf.py')
+md = md.extrude(3,1.)
+md.materials = matestar()
+md.materials.rheology_B = 3.15e8 * np.ones((md.mesh.numberofvertices,))
+md.materials.rheology_Ec = np.ones((md.mesh.numberofvertices,))
+md.materials.rheology_Es = 3. * np.ones((md.mesh.numberofvertices,))
+
+md = setflowequation(md,'FS','all')
+md.initialization.waterfraction = np.zeros((md.mesh.numberofvertices,))
+md.initialization.watercolumn = np.zeros((md.mesh.numberofvertices,))
+md.transient.isstressbalance = 0
+md.transient.ismasstransport = 0
+md.transient.issmb = 1
+md.transient.isthermal = 1
+md.transient.isgroundingline = 0
+md.thermal.isenthalpy = 1
+md.thermal.isdynamicbasalspc = 1
+md.cluster = generic('name',gethostname(),'np',3)
+md = solve(md,'Transient')
+
+#Fields and tolerances to track changes
+field_names = [
+       'Enthalpy1','Waterfraction1','Temperature1',
+       'Enthalpy2','Waterfraction2','Temperature2',
+       'Enthalpy3','Waterfraction3','Temperature3']
+field_tolerances = [
+       1e-12,1e-11,1e-12,
+       1e-12,1e-10,1e-12,
+       1e-12,1e-9,1e-12]
+field_values = [
+          md.results.TransientSolution[0].Enthalpy,
+          md.results.TransientSolution[0].Waterfraction,
+          md.results.TransientSolution[0].Temperature,
+          md.results.TransientSolution[1].Enthalpy,
+          md.results.TransientSolution[1].Waterfraction,
+          md.results.TransientSolution[1].Temperature,
+          md.results.TransientSolution[2].Enthalpy,
+          md.results.TransientSolution[2].Waterfraction,
+          md.results.TransientSolution[2].Temperature,
+          ]
Index: ../trunk-jpl/test/NightlyRun/test435.py
===================================================================
--- ../trunk-jpl/test/NightlyRun/test435.py     (nonexistent)
+++ ../trunk-jpl/test/NightlyRun/test435.py     (revision 22267)
@@ -0,0 +1,97 @@
+#Test Name: MISMIP3DHO
+import numpy as np
+from model import *
+from socket import gethostname
+from triangle import *
+from setmask import *
+from parameterize import *
+from setflowequation import *
+from solve import *
+
+md = triangle(model(),'../Exp/Square.exp',100000.)
+md = setmask(md,'../Exp/SquareShelf.exp','')
+md = parameterize(md,'../Par/SquareSheetShelf.py')
+md.initialization.vx[:] = 1.
+md.initialization.vy[:] = 1.
+md.geometry.thickness[:] = 500. - md.mesh.x / 10000.
+md.geometry.bed = -100. - md.mesh.x / 1000.
+md.geometry.base = -md.geometry.thickness * md.materials.rho_ice / md.materials.rho_water
+md.mask.groundedice_levelset = md.geometry.thickness + md.materials.rho_water / md.materials.rho_ice * md.geometry.bed
+pos = np.where(md.mask.groundedice_levelset >= 0)
+md.geometry.base[pos] = md.geometry.bed[pos]
+md.geometry.surface = md.geometry.base + md.geometry.thickness
+md = md.extrude(4,1.)
+md = setflowequation(md,'HO','all')
+
+#Boundary conditions:
+md.mask.ice_levelset = -np.ones((md.mesh.numberofvertices,))
+md.mask.ice_levelset[np.where(md.mesh.x == max(md.mesh.x))] = 0.
+md.stressbalance.spcvx[:] = float('Nan')
+md.stressbalance.spcvy[:] = float('Nan')
+md.stressbalance.spcvz[:] = float('Nan')
+posA = np.intersect1d(np.array(np.where(md.mesh.y < 1000000.1)),np.array(np.where(md.mesh.y > 999999.9)))
+posB = np.intersect1d(np.array(np.where(md.mesh.y < 0.1)),np.array(np.where(md.mesh.y > -0.1)))
+pos = np.unique(np.concatenate((posA,posB)))
+md.stressbalance.spcvy[pos] = 0.
+pos2 = np.intersect1d(np.array(np.where(md.mesh.x < 0.1)), np.array(np.where(md.mesh.x > -0.1)))
+md.stressbalance.spcvx[pos2] = 0.
+md.stressbalance.spcvy[pos2] = 0.
+
+md.materials.rheology_B = 1. / ((10**-25)**(1./3.)) * np.ones((md.mesh.numberofvertices,))
+md.materials.rheology_law = 'None'
+md.friction.coefficient[:] = np.sqrt(1e7) * np.ones((md.mesh.numberofvertices,))
+md.friction.p = 3. * np.ones((md.mesh.numberofelements,))
+md.smb.mass_balance[:] = 1.
+md.basalforcings.groundedice_melting_rate[:] = 0.
+md.basalforcings.floatingice_melting_rate[:] = 30.
+md.transient.isthermal = 0
+md.transient.isstressbalance = 1
+md.transient.isgroundingline = 1
+md.transient.ismasstransport = 1
+md.transient.issmb = 1
+md.transient.requested_outputs = ['default','BasalforcingsFloatingiceMeltingRate']
+md.groundingline.migration = 'SubelementMigration'
+md.timestepping.final_time = 30
+md.timestepping.time_step = 10
+
+md.cluster = generic('name',gethostname(),'np',3)
+md = solve(md,'Transient')
+
+#Fields and tolerances to track changes
+field_names     = [
+       'Bed1','Surface1','Thickness1','Floatingice1','Vx1','Vy1','Vz1','Pressure1','FloatingiceMeltingrate1',
+       'Bed2','Surface2','Thickness2','Floatingice2','Vx2','Vy2','Vz2','Pressure2','FloatingiceMeltingrate2',
+       'Bed3','Surface3','Thickness3','Floatingice3','Vx3','Vy3','Vz3','Pressure3','FloatingiceMeltingrate3']
+field_tolerances = [
+       2e-11,5e-12,2e-11,1e-11,5e-10,3e-08,6e-10,1e-13,1e-13,
+       3e-11,3e-11,9e-10,7e-11,6e-09,1e-07,1e-09,1e-10,1e-13,
+       1e-9,2e-08,7e-09,2e-7 ,1e-03,8e-04,2e-09,1e-10,1e-13]
+field_values = [
+       md.results.TransientSolution[0].Base,
+       md.results.TransientSolution[0].Surface,
+       md.results.TransientSolution[0].Thickness,
+       md.results.TransientSolution[0].MaskGroundediceLevelset,
+       md.results.TransientSolution[0].Vx,
+       md.results.TransientSolution[0].Vy,
+       md.results.TransientSolution[0].Vz,
+       md.results.TransientSolution[0].Pressure,
+       md.results.TransientSolution[0].BasalforcingsFloatingiceMeltingRate,
+       md.results.TransientSolution[1].Base,
+       md.results.TransientSolution[1].Surface,
+       md.results.TransientSolution[1].Thickness,
+       md.results.TransientSolution[1].MaskGroundediceLevelset,
+       md.results.TransientSolution[1].Vx,
+       md.results.TransientSolution[1].Vy,
+       md.results.TransientSolution[1].Vz,
+       md.results.TransientSolution[1].Pressure,
+       md.results.TransientSolution[1].BasalforcingsFloatingiceMeltingRate,
+       md.results.TransientSolution[2].Base,
+       md.results.TransientSolution[2].Surface,
+       md.results.TransientSolution[2].Thickness,
+       md.results.TransientSolution[2].MaskGroundediceLevelset,
+       md.results.TransientSolution[2].Vx,
+       md.results.TransientSolution[2].Vy,
+       md.results.TransientSolution[2].Vz,
+       md.results.TransientSolution[2].Pressure,
+       md.results.TransientSolution[2].BasalforcingsFloatingiceMeltingRate,
+       ]
Index: ../trunk-jpl/test/NightlyRun/test462.py
===================================================================
--- ../trunk-jpl/test/NightlyRun/test462.py     (nonexistent)
+++ ../trunk-jpl/test/NightlyRun/test462.py     (revision 22267)
@@ -0,0 +1,49 @@
+#Test Name: SquareSheetShelfAmrBamgField
+import numpy as np
+from model import *
+from socket import gethostname
+from triangle import *
+from setmask import *
+from parameterize import *
+from setflowequation import *
+from solve import *
+
+md = triangle(model(),'../Exp/Square.exp',150000.)
+md = setmask(md,'../Exp/SquareShelf.exp','')
+md = parameterize(md,'../Par/SquareSheetShelf.py')
+md = setflowequation(md,'SSA','all')
+md.cluster = generic('name',gethostname(),'np',3)
+md.transient.isstressbalance = 1
+md.transient.ismasstransport = 1
+md.transient.issmb = 0
+md.transient.isthermal = 0
+md.transient.isgroundingline = 0
+#amr bamg settings, just field
+md.amr.hmin = 10000
+md.amr.hmax = 100000
+md.amr.fieldname = 'Vel'
+md.amr.keepmetric = 1
+md.amr.gradation = 1.2
+md.amr.groundingline_resolution = 2000
+md.amr.groundingline_distance = 0
+md.amr.icefront_resolution = 1000
+md.amr.icefront_distance = 0
+md.amr.thicknesserror_resolution = 1000
+md.amr.thicknesserror_threshold = 0
+md.amr.deviatoricerror_resolution = 1000
+md.amr.deviatoricerror_threshold = 0
+md.transient.amr_frequency = 1
+md.timestepping.start_time = 0
+md.timestepping.final_time = 3
+md.timestepping.time_step = 1
+md = solve(md,'Transient')
+
+#Fields and tolerances to track changes
+field_names     = ['Vx','Vy','Vel','Pressure']
+field_tolerances = [1e-13,1e-13,1e-13,1e-13]
+field_values = [
+       md.results.TransientSolution[2].Vx,
+       md.results.TransientSolution[2].Vy,
+       md.results.TransientSolution[2].Vel,
+       md.results.TransientSolution[2].Pressure,
+       ]
Index: ../trunk-jpl/test/NightlyRun/test436.py
===================================================================
--- ../trunk-jpl/test/NightlyRun/test436.py     (nonexistent)
+++ ../trunk-jpl/test/NightlyRun/test436.py     (revision 22267)
@@ -0,0 +1,45 @@
+#Test Name: SquareSheetShelfSteaEnthalpyRheologiesHO
+import numpy as np
+from model import *
+from socket import gethostname
+from triangle import *
+from setmask import *
+from parameterize import *
+from setflowequation import *
+from solve import *
+
+md = triangle(model(),'../Exp/Square.exp',150000.)
+md = setmask(md,'../Exp/SquareShelf.exp','')
+md = parameterize(md,'../Par/SquareSheetShelf.py')
+md = md.extrude(3,2.)
+md = setflowequation(md,'HO','all')
+md.cluster = generic('name',gethostname(),'np',3)
+md.timestepping.time_step = 0.
+md.thermal.isenthalpy = 1
+md.initialization.waterfraction = np.zeros((md.mesh.numberofvertices,))
+md.initialization.watercolumn = np.zeros((md.mesh.numberofvertices,))
+
+#Go solve
+field_names = []
+field_tolerances = []
+field_values = []
+for i in ['LliboutryDuval', 'CuffeyTemperate']:
+       print ' '
+       print '====== Testing rheology law: ' + i + ' ====='
+
+       md.materials.rheology_law = i
+       md = solve(md,'Steadystate')
+       field_names += ['Vx'+i,'Vy'+i,'Vz'+i,'Vel'+i,'Pressure'+i,
+               'Temperature'+i,'Waterfraction'+i,'Enthalpy'+i]
+       field_tolerances += [2e-09,1e-09,1e-09,1e-09,1e-13,1e-10,6e-10,5e-10]
+       field_values += [
+               md.results.SteadystateSolution.Vx,
+               md.results.SteadystateSolution.Vy,
+               md.results.SteadystateSolution.Vz,
+               md.results.SteadystateSolution.Vel,
+               md.results.SteadystateSolution.Pressure,
+               md.results.SteadystateSolution.Temperature,
+               md.results.SteadystateSolution.Waterfraction,
+               md.results.SteadystateSolution.Enthalpy,
+               ]
+
Index: ../trunk-jpl/test/NightlyRun/test463.py
===================================================================
--- ../trunk-jpl/test/NightlyRun/test463.py     (nonexistent)
+++ ../trunk-jpl/test/NightlyRun/test463.py     (revision 22267)
@@ -0,0 +1,49 @@
+#Test Name: SquareSheetShelfAmrBamgGroundingline
+import numpy as np
+from model import *
+from socket import gethostname
+from triangle import *
+from setmask import *
+from parameterize import *
+from setflowequation import *
+from solve import *
+
+md = triangle(model(),'../Exp/Square.exp',150000.)
+md = setmask(md,'../Exp/SquareShelf.exp','')
+md = parameterize(md,'../Par/SquareSheetShelf.py')
+md = setflowequation(md,'SSA','all')
+md.cluster = generic('name',gethostname(),'np',3)
+md.transient.isstressbalance = 1
+md.transient.ismasstransport = 1
+md.transient.issmb = 0
+md.transient.isthermal = 0
+md.transient.isgroundingline = 0
+#amr bamg settings, just grounding line
+md.amr.hmin = 10000
+md.amr.hmax = 100000
+md.amr.fieldname = 'None'
+md.amr.keepmetric = 0
+md.amr.gradation = 1.2
+md.amr.groundingline_resolution = 12000
+md.amr.groundingline_distance = 100000
+md.amr.icefront_resolution = 1000
+md.amr.icefront_distance = 0
+md.amr.thicknesserror_resolution = 1000
+md.amr.thicknesserror_threshold = 0
+md.amr.deviatoricerror_resolution = 1000
+md.amr.deviatoricerror_threshold = 0
+md.transient.amr_frequency = 1
+md.timestepping.start_time = 0
+md.timestepping.final_time = 3
+md.timestepping.time_step = 1
+md = solve(md,'Transient')
+
+#Fields and tolerances to track changes
+field_names     = ['Vx','Vy','Vel','Pressure']
+field_tolerances = [1e-8,1e-8,1e-8,1e-8]
+field_values = [
+       md.results.TransientSolution[2].Vx,
+       md.results.TransientSolution[2].Vy,
+       md.results.TransientSolution[2].Vel,
+       md.results.TransientSolution[2].Pressure,
+       ]
Index: ../trunk-jpl/test/NightlyRun/test437.py
===================================================================
--- ../trunk-jpl/test/NightlyRun/test437.py     (nonexistent)
+++ ../trunk-jpl/test/NightlyRun/test437.py     (revision 22267)
@@ -0,0 +1,96 @@
+#Test Name: ThermalEnthBasalcondsTrans
+import numpy as np
+from model import *
+from socket import gethostname
+from triangle import *
+from setmask import *
+from parameterize import *
+from setflowequation import *
+from solve import *
+
+md = triangle(model(),'../Exp/Square.exp',300000.)
+md = setmask(md,'','')
+md = parameterize(md,'../Par/SquareThermal.py')
+
+h = 100.
+md.geometry.thickness = h * np.ones((md.mesh.numberofvertices,))
+md.geometry.base = -h * np.ones((md.mesh.numberofvertices,))
+md.geometry.surface = md.geometry.base + md.geometry.thickness
+
+md.extrude(41,2.)
+md = setflowequation(md,'HO','all')
+md.thermal.isenthalpy = True
+md.thermal.isdynamicbasalspc = True
+
+#Basal forcing
+Ts = 273.15-3.
+Tb = 273.15-1.
+Tsw = Tb
+qgeo = md.materials.thermalconductivity / max(md.geometry.thickness) * (Tb - Ts) #qgeo=kappa*(Tb-Ts)/H
+md.basalforcings.geothermalflux[np.where(md.mesh.vertexonbase)] = qgeo
+md.initialization.temperature = qgeo / md.materials.thermalconductivity * (md.geometry.surface - md.mesh.z) + Ts
+
+#Surface forcing
+pos = np.where(md.mesh.vertexonsurface)
+SPC_cold = float('NaN') * np.ones((md.mesh.numberofvertices,))
+SPC_warm = float('NaN') * np.ones((md.mesh.numberofvertices,))
+SPC_cold[pos] = Ts
+SPC_warm[pos] = Tsw
+md.thermal.spctemperature = SPC_cold
+md.timestepping.time_step = 5.
+t0 = 0.
+t1 = 10.
+t2 = 100.
+md.timestepping.final_time = 400.
+md.thermal.spctemperature = np.array([SPC_cold,SPC_cold,SPC_warm,SPC_warm,SPC_cold]).T
+md.thermal.spctemperature = np.vstack((md.thermal.spctemperature,[t0, t1-1, t1, t2-1, t2]))
+#print np.shape(md.thermal.spctemperature)
+#print md.thermal.spctemperature
+
+#Additional settings
+md.transient.ismasstransport = False
+md.transient.isstressbalance = False
+md.transient.issmb = True
+md.transient.isthermal = True
+md.thermal.stabilization = False
+
+#Go solve
+md.verbose = verbose(0)
+md.cluster = generic('name',gethostname(),'np',3)
+md = solve(md,'Transient')
+
+#Fields and tolerances to track changes
+field_names = ['Enthalpy1','Temperature1','Waterfraction1','BasalMeltingRate1','Watercolumn1',
+                  'Enthalpy2','Temperature2','Waterfraction2','BasalMeltingRate2','Watercolumn2',
+                  'Enthalpy3','Temperature3','Waterfraction3','BasalMeltingRate3','Watercolumn3',
+                  'Enthalpy4','Temperature4','Waterfraction4','BasalMeltingRate4','Watercolumn4']
+field_tolerances = [1.e-10,1.e-10,1.e-10,1.e-9,1.e-10,
+                   1.e-10,1.e-10,1.e-10,2.e-9,2.e-10,
+                   1.e-10,1.e-10,1.e-10,2.e-9,1.e-10,
+                   1.e-10,1.e-10,1.e-10,2.e-9,1.e-10]
+i1 = 0
+i2 = int(np.ceil(t2 / md.timestepping.time_step) + 2)-1
+i3 = int(np.ceil(md.timestepping.final_time / (2. * md.timestepping.time_step)))-1
+i4 = np.shape(md.results.TransientSolution)[0]-1
+field_values = [
+       md.results.TransientSolution[i1].Enthalpy,
+       md.results.TransientSolution[i1].Temperature,
+       md.results.TransientSolution[i1].Waterfraction,
+       md.results.TransientSolution[i1].BasalforcingsGroundediceMeltingRate,
+       md.results.TransientSolution[i1].Watercolumn,
+       md.results.TransientSolution[i2].Enthalpy,
+       md.results.TransientSolution[i2].Temperature,
+       md.results.TransientSolution[i2].Waterfraction,
+       md.results.TransientSolution[i2].BasalforcingsGroundediceMeltingRate,
+       md.results.TransientSolution[i2].Watercolumn,
+       md.results.TransientSolution[i3].Enthalpy,
+       md.results.TransientSolution[i3].Temperature,
+       md.results.TransientSolution[i3].Waterfraction,
+       md.results.TransientSolution[i3].BasalforcingsGroundediceMeltingRate,
+       md.results.TransientSolution[i3].Watercolumn,
+       md.results.TransientSolution[i4].Enthalpy,
+       md.results.TransientSolution[i4].Temperature,
+       md.results.TransientSolution[i4].Waterfraction,
+       md.results.TransientSolution[i4].BasalforcingsGroundediceMeltingRate,
+       md.results.TransientSolution[i4].Watercolumn,
+       ]
Index: ../trunk-jpl/test/NightlyRun/test464.py
===================================================================
--- ../trunk-jpl/test/NightlyRun/test464.py     (nonexistent)
+++ ../trunk-jpl/test/NightlyRun/test464.py     (revision 22267)
@@ -0,0 +1,49 @@
+#Test Name: SquareSheetShelfAmrBamgIceFront
+import numpy as np
+from model import *
+from socket import gethostname
+from triangle import *
+from setmask import *
+from parameterize import *
+from setflowequation import *
+from solve import *
+
+md = triangle(model(),'../Exp/Square.exp',150000.)
+md = setmask(md,'../Exp/SquareShelf.exp','')
+md = parameterize(md,'../Par/SquareSheetShelf.py')
+md = setflowequation(md,'SSA','all')
+md.cluster = generic('name',gethostname(),'np',3)
+md.transient.isstressbalance = 1
+md.transient.ismasstransport = 1
+md.transient.issmb = 0
+md.transient.isthermal = 0
+md.transient.isgroundingline = 0
+#amr bamg settings, just ice front
+md.amr.hmin = 10000
+md.amr.hmax = 100000
+md.amr.fieldname = 'None'
+md.amr.keepmetric = 0
+md.amr.gradation = 1.2
+md.amr.groundingline_resolution = 12000
+md.amr.groundingline_distance = 0
+md.amr.icefront_resolution = 12000
+md.amr.icefront_distance = 100000
+md.amr.thicknesserror_resolution = 1000
+md.amr.thicknesserror_threshold = 0
+md.amr.deviatoricerror_resolution = 1000
+md.amr.deviatoricerror_threshold = 0
+md.transient.amr_frequency = 1
+md.timestepping.start_time = 0
+md.timestepping.final_time = 3
+md.timestepping.time_step = 1
+md = solve(md,'Transient')
+
+#Fields and tolerances to track changes
+field_names     = ['Vx','Vy','Vel','Pressure']
+field_tolerances = [1e-13,1e-13,1e-13,1e-13]
+field_values = [
+       md.results.TransientSolution[2].Vx,
+       md.results.TransientSolution[2].Vy,
+       md.results.TransientSolution[2].Vel,
+       md.results.TransientSolution[2].Pressure,
+       ]
Index: ../trunk-jpl/test/NightlyRun/test438.py
===================================================================
--- ../trunk-jpl/test/NightlyRun/test438.py     (nonexistent)
+++ ../trunk-jpl/test/NightlyRun/test438.py     (revision 22267)
@@ -0,0 +1,53 @@
+#Test Name: TransientFrictionWaterlayer2D
+import numpy as np
+from model import *
+from socket import gethostname
+from triangle import *
+from setmask import *
+from parameterize import *
+from setflowequation import *
+from solve import *
+from frictionwaterlayer import *
+
+md = triangle(model(),'../Exp/Square.exp',150000.)
+md = setmask(md,'../Exp/SquareShelf.exp','')
+md = parameterize(md,'../Par/SquareSheetShelf.py')
+md = setflowequation(md,'SSA','all')
+md.friction = frictionwaterlayer(md)
+md.friction.water_layer = np.zeros((md.mesh.numberofvertices+1,2))
+md.friction.water_layer[:,1] = 1.
+md.friction.water_layer[md.mesh.numberofvertices,:] = [1.,2.]
+md.friction.f = 0.8
+md.cluster = generic('name',gethostname(),'np',3)
+md = solve(md,'Transient')
+
+#Fields and tolerances to track changes
+field_names =['Vx1','Vy1','Vel1','Pressure1','Bed1','Surface1','Thickness1',
+             'Vx2','Vy2','Vel2','Pressure2','Bed2','Surface2','Thickness2',
+             'Vx3','Vy3','Vel3','Pressure3','Bed3','Surface3','Thickness3']
+field_tolerances=[1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,
+                 1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,
+                 1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13]
+field_values = [
+       md.results.TransientSolution[0].Vx,
+       md.results.TransientSolution[0].Vy,
+       md.results.TransientSolution[0].Vel,
+       md.results.TransientSolution[0].Pressure,
+       md.results.TransientSolution[0].Base,
+       md.results.TransientSolution[0].Surface,
+       md.results.TransientSolution[0].Thickness,
+       md.results.TransientSolution[1].Vx,
+       md.results.TransientSolution[1].Vy,
+       md.results.TransientSolution[1].Vel,
+       md.results.TransientSolution[1].Pressure,
+       md.results.TransientSolution[1].Base,
+       md.results.TransientSolution[1].Surface,
+       md.results.TransientSolution[1].Thickness,
+       md.results.TransientSolution[2].Vx,
+       md.results.TransientSolution[2].Vy,
+       md.results.TransientSolution[2].Vel,
+       md.results.TransientSolution[2].Pressure,
+       md.results.TransientSolution[2].Base,
+       md.results.TransientSolution[2].Surface,
+       md.results.TransientSolution[2].Thickness
+       ]
Index: ../trunk-jpl/test/NightlyRun/test465.py
===================================================================
--- ../trunk-jpl/test/NightlyRun/test465.py     (nonexistent)
+++ ../trunk-jpl/test/NightlyRun/test465.py     (revision 22267)
@@ -0,0 +1,49 @@
+#Test Name: SquareSheetShelfAmrBamgAll
+import numpy as np
+from model import *
+from socket import gethostname
+from triangle import *
+from setmask import *
+from parameterize import *
+from setflowequation import *
+from solve import *
+
+md = triangle(model(),'../Exp/Square.exp',150000.)
+md = setmask(md,'../Exp/SquareShelf.exp','')
+md = parameterize(md,'../Par/SquareSheetShelf.py')
+md = setflowequation(md,'SSA','all')
+md.cluster = generic('name',gethostname(),'np',3)
+md.transient.isstressbalance = 1
+md.transient.ismasstransport = 1
+md.transient.issmb = 0
+md.transient.isthermal = 0
+md.transient.isgroundingline = 0
+#amr bamg settings, field, grounding line and ice front
+md.amr.hmin = 20000
+md.amr.hmax = 100000
+md.amr.fieldname = 'Vel'
+md.amr.keepmetric = 0
+md.amr.gradation = 1.2
+md.amr.groundingline_resolution = 10000
+md.amr.groundingline_distance = 100000
+md.amr.icefront_resolution = 15000
+md.amr.icefront_distance = 100000
+md.amr.thicknesserror_resolution = 1000
+md.amr.thicknesserror_threshold = 0
+md.amr.deviatoricerror_resolution = 1000
+md.amr.deviatoricerror_threshold = 0
+md.transient.amr_frequency = 1
+md.timestepping.start_time = 0
+md.timestepping.final_time = 3
+md.timestepping.time_step = 1
+md = solve(md,'Transient')
+
+#Fields and tolerances to track changes
+field_names     = ['Vx','Vy','Vel','Pressure']
+field_tolerances = [1e-8,1e-8,1e-8,1e-8]
+field_values = [
+       md.results.TransientSolution[2].Vx,
+       md.results.TransientSolution[2].Vy,
+       md.results.TransientSolution[2].Vel,
+       md.results.TransientSolution[2].Pressure,
+       ]
Index: ../trunk-jpl/test/NightlyRun/test439.py
===================================================================
--- ../trunk-jpl/test/NightlyRun/test439.py     (nonexistent)
+++ ../trunk-jpl/test/NightlyRun/test439.py     (revision 22267)
@@ -0,0 +1,54 @@
+#Test Name: TransientFrictionWaterlayer3D
+import numpy as np
+from model import *
+from socket import gethostname
+from triangle import *
+from setmask import *
+from parameterize import *
+from setflowequation import *
+from solve import *
+from frictionwaterlayer import *
+
+md = triangle(model(),'../Exp/Square.exp',150000.)
+md = setmask(md,'../Exp/SquareShelf.exp','')
+md = parameterize(md,'../Par/SquareSheetShelf.py')
+md = md.extrude(4,1.)
+md = setflowequation(md,'HO','all')
+md.friction = frictionwaterlayer(md)
+md.friction.water_layer = np.zeros((md.mesh.numberofvertices+1,2))
+md.friction.water_layer[:,1] = 1.
+md.friction.water_layer[md.mesh.numberofvertices,:] = [1.,2.]
+md.friction.f = 0.8
+md.cluster = generic('name',gethostname(),'np',3)
+md = solve(md,'Transient')
+
+#Fields and tolerances to track changes
+field_names =['Vx1','Vy1','Vel1','Pressure1','Bed1','Surface1','Thickness1',
+             'Vx2','Vy2','Vel2','Pressure2','Bed2','Surface2','Thickness2',
+             'Vx3','Vy3','Vel3','Pressure3','Bed3','Surface3','Thickness3']
+field_tolerances = [1e-09,1e-09,1e-09,1e-09,1e-09,1e-09,1e-09,
+                   1e-09,1e-09,1e-09,1e-09,1e-09,1e-09,1e-09,
+                   1e-09,1e-09,1e-09,1e-09,1e-09,1e-09,1e-09]
+field_values = [
+       md.results.TransientSolution[0].Vx,
+       md.results.TransientSolution[0].Vy,
+       md.results.TransientSolution[0].Vel,
+       md.results.TransientSolution[0].Pressure,
+       md.results.TransientSolution[0].Base,
+       md.results.TransientSolution[0].Surface,
+       md.results.TransientSolution[0].Thickness,
+       md.results.TransientSolution[1].Vx,
+       md.results.TransientSolution[1].Vy,
+       md.results.TransientSolution[1].Vel,
+       md.results.TransientSolution[1].Pressure,
+       md.results.TransientSolution[1].Base,
+       md.results.TransientSolution[1].Surface,
+       md.results.TransientSolution[1].Thickness,
+       md.results.TransientSolution[2].Vx,
+       md.results.TransientSolution[2].Vy,
+       md.results.TransientSolution[2].Vel,
+       md.results.TransientSolution[2].Pressure,
+       md.results.TransientSolution[2].Base,
+       md.results.TransientSolution[2].Surface,
+       md.results.TransientSolution[2].Thickness
+       ]
Index: ../trunk-jpl/test/NightlyRun/test340.py
===================================================================
--- ../trunk-jpl/test/NightlyRun/test340.py     (nonexistent)
+++ ../trunk-jpl/test/NightlyRun/test340.py     (revision 22267)
@@ -0,0 +1,46 @@
+#Test Name: SquareSheetConstrainedSmbGradientsEla3d
+import numpy as np
+from model import *
+from socket import gethostname
+from triangle import *
+from setmask import *
+from parameterize import *
+from setflowequation import *
+from solve import *
+from taoinversion import *
+
+md = triangle(model(),'../Exp/Square.exp',200000.)
+md = setmask(md,'','')
+md = parameterize(md,'../Par/SquareSheetConstrained.py')
+md.extrude(3,1.)
+md = setflowequation(md,'HO','all')
+
+#control parameters
+md.inversion = taoinversion()
+md.inversion.iscontrol = 1
+md.inversion.control_parameters = ['FrictionCoefficient']
+md.inversion.min_parameters = 1. * np.ones((md.mesh.numberofvertices,))
+md.inversion.max_parameters = 200. * np.ones((md.mesh.numberofvertices,))
+md.inversion.maxsteps = 2
+md.inversion.maxiter = 6
+md.inversion.cost_functions = [102,501]
+md.inversion.cost_functions_coefficients = np.ones((md.mesh.numberofvertices,2))
+md.inversion.cost_functions_coefficients[:,1] = 2. * 1e-7
+md.inversion.vx_obs = md.initialization.vx
+md.inversion.vy_obs = md.initialization.vy
+
+md.cluster = generic('name',gethostname(),'np',3)
+md = solve(md,'Stressbalance')
+
+#Fields and tolerances to track changes
+field_names      = ['Gradient','Misfits','FrictionCoefficient','Pressure','Vel','Vx','Vy']
+field_tolerances = [3e-08,1e-07,5e-10,1e-10,1e-09,1e-09,1e-09]
+field_values = [
+       md.results.StressbalanceSolution.Gradient1,
+       md.results.StressbalanceSolution.J,
+       md.results.StressbalanceSolution.FrictionCoefficient,
+       md.results.StressbalanceSolution.Pressure,
+       md.results.StressbalanceSolution.Vel,
+       md.results.StressbalanceSolution.Vx,
+       md.results.StressbalanceSolution.Vy
+]
Index: ../trunk-jpl/test/NightlyRun/test430.py
===================================================================
--- ../trunk-jpl/test/NightlyRun/test430.py     (nonexistent)
+++ ../trunk-jpl/test/NightlyRun/test430.py     (revision 22267)
@@ -0,0 +1,95 @@
+#Test Name: MISMIP3D
+import numpy as np
+from model import *
+from socket import gethostname
+from triangle import *
+from setmask import *
+from parameterize import *
+from setflowequation import *
+from solve import *
+
+md = triangle(model(),'../Exp/Square.exp',100000.)
+md = setmask(md,'../Exp/SquareShelf.exp','')
+md = parameterize(md,'../Par/SquareSheetShelf.py')
+md.initialization.vx[:] = 1.
+md.initialization.vy[:] = 1.
+md.geometry.thickness[:] = 500. - md.mesh.x / 10000.
+md.geometry.bed = -100. - md.mesh.x / 1000.
+md.geometry.base = -md.geometry.thickness * md.materials.rho_ice / md.materials.rho_water
+md.mask.groundedice_levelset = md.geometry.thickness + md.materials.rho_water / md.materials.rho_ice * md.geometry.bed
+pos = np.where(md.mask.groundedice_levelset >= 0.)
+md.geometry.base[pos] = md.geometry.bed[pos]
+md.geometry.surface = md.geometry.base + md.geometry.thickness
+md = setflowequation(md,'SSA','all')
+
+#Boundary conditions:
+md.mask.ice_levelset = -np.ones((md.mesh.numberofvertices,))
+md.mask.ice_levelset[np.where(md.mesh.x == max(md.mesh.x))] = 0.
+md.stressbalance.spcvx[:] = float('NaN')
+md.stressbalance.spcvy[:] = float('NaN')
+md.stressbalance.spcvz[:] = float('NaN')
+posA = np.intersect1d(np.array(np.where(md.mesh.y < 1000000.1)),np.array(np.where(md.mesh.y > 999999.9)))
+posB = np.intersect1d(np.array(np.where(md.mesh.y < 0.1)),np.array(np.where(md.mesh.y > -0.1)))
+pos = np.unique(np.concatenate((posA,posB)))
+md.stressbalance.spcvy[pos] = 0.
+pos2 = np.intersect1d(np.array(np.where(md.mesh.x < 0.1)), np.array(np.where(md.mesh.x > -0.1)))
+md.stressbalance.spcvx[pos2] = 0.
+md.stressbalance.spcvy[pos2] = 0.
+
+md.materials.rheology_B = 1. / ((10**-25)**(1. / 3.)) * np.ones((md.mesh.numberofvertices,))
+md.materials.rheology_law = 'None'
+md.friction.coefficient[:] = np.sqrt(10**7) * np.ones((md.mesh.numberofvertices,))
+md.friction.p = 3. * np.ones((md.mesh.numberofelements,))
+md.smb.mass_balance[:] = 1.
+md.basalforcings.groundedice_melting_rate[:] = 0.
+md.basalforcings.floatingice_melting_rate[:] = 30.
+md.transient.isthermal = 0
+md.transient.isstressbalance = 1
+md.transient.isgroundingline = 1
+md.transient.ismasstransport = 1
+md.transient.issmb = 1
+md.transient.requested_outputs = ['default','BasalforcingsFloatingiceMeltingRate']
+md.groundingline.migration = 'SubelementMigration'
+md.timestepping.final_time = 30
+md.timestepping.time_step = 10
+
+md.cluster = generic('name',gethostname(),'np',3)
+md = solve(md,'Transient')
+#print md.results.TransientSolution[0].BasalforcingsFloatingiceMeltingRate
+#print md.results.TransientSolution[1].BasalforcingsFloatingiceMeltingRate
+#print md.results.TransientSolution[2].BasalforcingsFloatingiceMeltingRate
+
+#Fields and tolerances to track changes
+field_names     = [
+       'Bed1','Surface1','Thickness1','Floatingice1','Vx1','Vy1','Pressure1','FloatingiceMeltingrate1',
+       'Bed2','Surface2','Thickness2','Floatingice2','Vx2','Vy2','Pressure2','FloatingiceMeltingrate2',
+       'Bed3','Surface3','Thickness3','Floatingice3','Vx3','Vy3','Pressure3','FloatingiceMeltingrate3']
+field_tolerances = [2e-11,5e-12,2e-11,1e-11,5e-10,1e-08,1e-13,1e-13,
+       3e-11,3e-11,9e-10,7e-11,1e-09,5e-08,1e-10,1e-13,
+       1e-10,3e-11,1e-10,7e-11,1e-09,5e-08,1e-10,1e-13]
+field_values = [
+       md.results.TransientSolution[0].Base,
+       md.results.TransientSolution[0].Surface,
+       md.results.TransientSolution[0].Thickness,
+       md.results.TransientSolution[0].MaskGroundediceLevelset,
+       md.results.TransientSolution[0].Vx,
+       md.results.TransientSolution[0].Vy,
+       md.results.TransientSolution[0].Pressure,
+       md.results.TransientSolution[0].BasalforcingsFloatingiceMeltingRate,
+       md.results.TransientSolution[1].Base,
+       md.results.TransientSolution[1].Surface,
+       md.results.TransientSolution[1].Thickness,
+       md.results.TransientSolution[1].MaskGroundediceLevelset,
+       md.results.TransientSolution[1].Vx,
+       md.results.TransientSolution[1].Vy,
+       md.results.TransientSolution[1].Pressure,
+       md.results.TransientSolution[1].BasalforcingsFloatingiceMeltingRate,
+       md.results.TransientSolution[2].Base,
+       md.results.TransientSolution[2].Surface,
+       md.results.TransientSolution[2].Thickness,
+       md.results.TransientSolution[2].MaskGroundediceLevelset,
+       md.results.TransientSolution[2].Vx,
+       md.results.TransientSolution[2].Vy,
+       md.results.TransientSolution[2].Pressure,
+       md.results.TransientSolution[2].BasalforcingsFloatingiceMeltingRate,
+       ]
Index: ../trunk-jpl/test/NightlyRun/test808.py
===================================================================
--- ../trunk-jpl/test/NightlyRun/test808.py     (nonexistent)
+++ ../trunk-jpl/test/NightlyRun/test808.py     (revision 22267)
@@ -0,0 +1,76 @@
+#Test Name: SquareShelfLevelsetCalvingSSA2dMinThickness
+from model import *
+from socket import gethostname
+from triangle import *
+from setmask import *
+from parameterize import *
+from setflowequation import *
+from solve import *
+import numpy as np
+from calvingminthickness import *
+
+md = triangle(model(),'../Exp/Square.exp',30000.)
+md = setmask(md,'all','')
+md = parameterize(md,'../Par/SquareShelf.py')
+md = setflowequation(md,'SSA','all')
+md.cluster = generic('name',gethostname(),'np',3)
+
+x = md.mesh.x
+xmin = min(x)
+xmax = max(x)
+Lx = (xmax - xmin)
+alpha = 2. / 3.
+md.mask.ice_levelset = -1 + 2 * (md.mesh.y > 9e5)
+
+md.timestepping.time_step = 1
+md.timestepping.final_time = 30
+
+#Transient
+md.transient.isstressbalance = 1
+md.transient.ismasstransport = 1
+md.transient.issmb = 1
+md.transient.isthermal = 0
+md.transient.isgroundingline = 0
+md.transient.isgia = 0
+md.transient.ismovingfront = 1
+
+md.calving = calvingminthickness()
+md.calving.min_thickness = 400
+md.calving.meltingrate = np.zeros((md.mesh.numberofvertices,))
+md.levelset.spclevelset = float('NaN')* np.ones((md.mesh.numberofvertices,))
+md.levelset.reinit_frequency = 1
+
+md = solve(md,'Transient')
+
+#Fields and tolerances to track changes
+field_names = [
+       'Vx1','Vy1','Vel1','Pressure1','Thickness1','Surface1','MaskIceLevelset1'
+       'Vx2','Vy2','Vel2','Pressure2','Thickness2','Surface2','MaskIceLevelset2'
+       'Vx3','Vy3','Vel3','Pressure3','Thickness3','Surface3','MaskIceLevelset3']
+field_tolerances = [
+       1e-8,1e-8,1e-8,1e-9,1e-9,1e-9,3e-9,
+       1e-8,1e-8,1e-8,1e-9,1e-9,1e-9,3e-9,
+       1e-8,1e-8,1e-8,1e-9,1e-9,1e-9,3e-9]
+field_values = [
+       md.results.TransientSolution[0].Vx,
+       md.results.TransientSolution[0].Vy,
+       md.results.TransientSolution[0].Vel,
+       md.results.TransientSolution[0].Pressure,
+       md.results.TransientSolution[0].Thickness,
+       md.results.TransientSolution[0].Surface,
+       md.results.TransientSolution[0].MaskIceLevelset,
+       md.results.TransientSolution[1].Vx,
+       md.results.TransientSolution[1].Vy,
+       md.results.TransientSolution[1].Vel,
+       md.results.TransientSolution[1].Pressure,
+       md.results.TransientSolution[1].Thickness,
+       md.results.TransientSolution[1].Surface,
+       md.results.TransientSolution[1].MaskIceLevelset,
+       md.results.TransientSolution[2].Vx,
+       md.results.TransientSolution[2].Vy,
+       md.results.TransientSolution[2].Vel,
+       md.results.TransientSolution[2].Pressure,
+       md.results.TransientSolution[2].Thickness,
+       md.results.TransientSolution[2].Surface,
+       md.results.TransientSolution[2].MaskIceLevelset
+       ]
Index: ../trunk-jpl/test/NightlyRun/test2424.py
===================================================================
--- ../trunk-jpl/test/NightlyRun/test2424.py    (nonexistent)
+++ ../trunk-jpl/test/NightlyRun/test2424.py    (revision 22267)
@@ -0,0 +1,51 @@
+#Test Name: SquareSheetShelfGroundingLine2dAggressive. From test424, with sea level increasing.
+import numpy as np
+from model import *
+from socket import gethostname
+from triangle import *
+from setmask import *
+from parameterize import *
+from setflowequation import *
+from solve import *
+from newforcing import *
+
+md = triangle(model(),'../Exp/Square.exp',150000.)
+md = setmask(md,'../Exp/SquareShelf.exp','')
+md = parameterize(md,'../Par/SquareSheetShelf.py')
+md = setflowequation(md,'SSA','all')
+md.initialization.vx[:] = 0.
+md.initialization.vy[:] = 0.
+md.smb.mass_balance[:] = 0.
+
+md.geometry.base = -700. - np.abs(md.mesh.y-500000.) / 1000.
+md.geometry.bed = -700. - np.abs(md.mesh.y-500000.) / 1000.
+md.geometry.thickness[:] = 1000.
+md.geometry.surface = md.geometry.base + md.geometry.thickness
+
+md.transient.isstressbalance = 0
+md.transient.isgroundingline = 1
+md.transient.isthermal = 0
+md.groundingline.migration = 'AggressiveMigration'
+md.transient.requested_outputs = ['IceVolume','IceVolumeAboveFloatation','Sealevel']
+
+md.timestepping.time_step = .1
+md.slr.sealevel = newforcing(md.timestepping.start_time, md.timestepping.final_time,
+                            md.timestepping.time_step, -200., 200., md.mesh.numberofvertices)
+
+md.cluster = generic('name',gethostname(),'np',3)
+md = solve(md,'Transient')
+
+#we are checking that the grounding line position is near the theorical one, which is the 0 contour level 
+#of surface - sealevel - (1-di)* thickness 
+
+nsteps = len(md.results.TransientSolution)
+field_names = []
+field_tolerances = []
+field_values = []
+#time is off by the year constant
+for i in range(nsteps):
+       field_names.append('Time-' + str(md.results.TransientSolution[i].time) + 
+               '-yr-ice_levelset-S-sl-(1-di)*H')
+       field_tolerances.append(1e-12)
+       field_values.append(md.results.TransientSolution[i].MaskGroundediceLevelset.reshape(-1,) - (md.geometry.surface - md.results.TransientSolution[i].Sealevel.reshape(-1,) - (1 - md.materials.rho_ice / md.materials.rho_water) * md.geometry.thickness))
+
Index: ../trunk-jpl/test/NightlyRun/test260.py
===================================================================
--- ../trunk-jpl/test/NightlyRun/test260.py     (nonexistent)
+++ ../trunk-jpl/test/NightlyRun/test260.py     (revision 22267)
@@ -0,0 +1,31 @@
+#Test Name: SquareShelfStressSSA2dEnhanced
+import numpy as np
+from model import *
+from socket import gethostname
+from triangle import *
+from setmask import *
+from parameterize import *
+from setflowequation import *
+from solve import *
+from matenhancedice import *
+
+md = triangle(model(),'../Exp/Square.exp',150000.)
+md = setmask(md,'all','')
+md.materials = matenhancedice()
+md.materials.rheology_B = 3.15e8 * np.ones(md.mesh.numberofvertices,)
+md.materials.rheology_n = 3 * np.ones(md.mesh.numberofelements,)
+md.materials.rheology_E = np.ones(md.mesh.numberofvertices,)
+md = parameterize(md,'../Par/SquareShelf.py')
+md = setflowequation(md,'SSA','all')
+md.cluster = generic('name',gethostname(),'np',3)
+md = solve(md,'Stressbalance')
+
+#Fields and tolerances to track changes
+field_names      = ['Vx','Vy','Vel','Pressure']
+field_tolerances = [1e-13,1e-13,1e-13,1e-13]
+field_values = [
+       md.results.StressbalanceSolution.Vx,
+       md.results.StressbalanceSolution.Vy,
+       md.results.StressbalanceSolution.Vel,
+       md.results.StressbalanceSolution.Pressure,
+       ]
Index: ../trunk-jpl/test/NightlyRun/test350.py
===================================================================
--- ../trunk-jpl/test/NightlyRun/test350.py     (nonexistent)
+++ ../trunk-jpl/test/NightlyRun/test350.py     (revision 22267)
@@ -0,0 +1,94 @@
+#Test Name: SquareSheetHydrologySommers
+from operator import itemgetter
+import numpy as np
+from model import *
+from socket import gethostname
+from triangle import *
+from setmask import *
+from parameterize import *
+from setflowequation import *
+from solve import *
+from frictionsommers import *
+from hydrologysommers import *
+from transient import *
+
+md = triangle(model(),'../Exp/Square.exp',50000.)
+md.mesh.x = md.mesh.x / 1000
+md.mesh.y = md.mesh.y / 1000
+md = setmask(md,'','')
+md = parameterize(md,'../Par/SquareSheetConstrained.py')
+md.transient = transient().deactivateall()
+md.transient.ishydrology = 1
+md = setflowequation(md,'SSA','all')
+md.cluster = generic('name',gethostname(),'np',2)
+
+#Use hydroogy coupled friciton law
+md.friction = frictionsommers(md)
+
+#Change hydrology class to Sommers' model
+md.hydrology = hydrologysommers()
+
+#Change geometry
+md.geometry.base = -.02 * md.mesh.x + 20.
+md.geometry.thickness = 300. * np.ones((md.mesh.numberofvertices,))
+md.geometry.bed = md.geometry.base
+md.geometry.surface = md.geometry.bed + md.geometry.thickness
+
+#define the initial water head as being such that the water pressure is 50% of the ice overburden pressure
+md.hydrology.head = 0.5 * md.materials.rho_ice / md.materials.rho_freshwater * md.geometry.thickness + md.geometry.base
+md.hydrology.gap_height = 0.01 * np.ones((md.mesh.numberofelements,))
+md.hydrology.bump_spacing = 2 * np.ones((md.mesh.numberofelements,))
+md.hydrology.bump_height = 0.05 * np.ones((md.mesh.numberofelements,))
+md.hydrology.englacial_input = 0.5 * np.ones((md.mesh.numberofvertices,))
+md.hydrology.reynolds= 1000. * np.ones((md.mesh.numberofelements,))
+md.hydrology.spchead = float('NaN') * np.ones((md.mesh.numberofvertices,))
+pos = np.intersect1d(np.array(np.where(md.mesh.vertexonboundary)), np.array(np.where(md.mesh.x == 1000)))
+md.hydrology.spchead[pos] = md.geometry.base[pos]
+
+#Define velocity
+md.initialization.vx = 1e-6 * md.constants.yts * np.ones((md.mesh.numberofvertices,))
+md.initialization.vy = np.zeros((md.mesh.numberofvertices,))
+
+md.timestepping.time_step = 3. * 3600. / md.constants.yts
+md.timestepping.final_time = .5 / 365.
+md.materials.rheology_B = (5e-25)**(-1. / 3.) * np.ones((md.mesh.numberofvertices,))
+
+#Add one moulin and Neumann BC, varying in time
+a = np.sqrt((md.mesh.x - 500.)**2 + (md.mesh.y - 500.)**2)
+pos = min(enumerate(a), key=itemgetter(1))[0]
+time = np.arange(0,md.timestepping.final_time+1,md.timestepping.time_step)
+md.hydrology.moulin_input = np.zeros((md.mesh.numberofvertices+1,np.size(time)))
+md.hydrology.moulin_input[-1,:] = time
+md.hydrology.moulin_input[pos,:] = 5. * (1. - np.sin(2. * np.pi / (1. / 365.) * time))
+md.hydrology.neumannflux = np.zeros((md.mesh.numberofelements+1,np.size(time)))
+md.hydrology.neumannflux[-1,:] = time
+segx = md.mesh.x[md.mesh.segments[:,0]-1]
+segy = md.mesh.y[md.mesh.segments[:,0]-1]
+posA = np.intersect1d(np.intersect1d(np.array(np.where(segx < 1.)),np.array(np.where(segy > 400.))), np.array(np.where(segy < 600.)))
+pos = (md.mesh.segments[posA]-1)[:,2]
+md.hydrology.neumannflux[pos,:] = np.tile(0.05*(1.-np.sin(2.*np.pi/(1./365.)*time)),(len(pos),1))
+
+md = solve(md,'Transient')
+
+#Fields and tolerances to track changes
+field_names = [
+       'HydrologyHead1','HydrologyGapHeight1',
+       'HydrologyHead2','HydrologyGapHeight2',
+       'HydrologyHead3','HydrologyGapHeight3',
+       'HydrologyHead4','HydrologyGapHeight4']
+field_tolerances = [
+       1e-13, 1e-13,
+       1e-13, 1e-13,
+       1e-13, 1e-13,
+       1e-13, 1e-12]
+field_values = [
+       md.results.TransientSolution[0].HydrologyHead,
+       md.results.TransientSolution[0].HydrologyGapHeight,
+       md.results.TransientSolution[1].HydrologyHead,
+       md.results.TransientSolution[1].HydrologyGapHeight,
+       md.results.TransientSolution[2].HydrologyHead,
+       md.results.TransientSolution[2].HydrologyGapHeight,
+       md.results.TransientSolution[3].HydrologyHead,
+       md.results.TransientSolution[3].HydrologyGapHeight
+       ]
+
Index: ../trunk-jpl/test/NightlyRun/test243.py
===================================================================
--- ../trunk-jpl/test/NightlyRun/test243.py     (nonexistent)
+++ ../trunk-jpl/test/NightlyRun/test243.py     (revision 22267)
@@ -0,0 +1,72 @@
+#Test Name: SquareShelfSMBGemb
+import numpy as np
+import scipy.io as spio
+from model import *
+from socket import gethostname
+from triangle import *
+from setmask import *
+from parameterize import *
+from setflowequation import *
+from solve import *
+from SMBgemb import *
+
+md = triangle(model(),'../Exp/Square.exp',200000.)
+md = setmask(md,'all','')
+md = parameterize(md,'../Par/SquareShelf.py')
+md = setflowequation(md,'SSA','all')
+md.materials.rho_ice = 910
+md.cluster = generic('name',gethostname(),'np',3)
+
+#Use of Gemb method for SMB computation
+md.smb = SMBgemb()
+md.smb.setdefaultparameters(md.mesh,md.geometry)
+
+#load hourly surface forcing date from 1979 to 2009:
+inputs = spio.loadmat('../Data/gemb_input.mat',squeeze_me=True)
+
+#setup the inputs: 
+md.smb.Ta = np.append(np.tile(np.conjugate(inputs['Ta0']),(md.mesh.numberofelements,1)),np.conjugate([inputs['dateN']]),axis=0)
+md.smb.V = np.append(np.tile(np.conjugate(inputs['V0']),(md.mesh.numberofelements,1)),np.conjugate([inputs['dateN']]),axis=0)
+md.smb.dswrf = np.append(np.tile(np.conjugate(inputs['dsw0']),(md.mesh.numberofelements,1)),np.conjugate([inputs['dateN']]),axis=0)
+md.smb.dlwrf = np.append(np.tile(np.conjugate(inputs['dlw0']),(md.mesh.numberofelements,1)),np.conjugate([inputs['dateN']]),axis=0)
+md.smb.P = np.append(np.tile(np.conjugate(inputs['P0']),(md.mesh.numberofelements,1)),np.conjugate([inputs['dateN']]),axis=0)
+md.smb.eAir = np.append(np.tile(np.conjugate(inputs['eAir0']),(md.mesh.numberofelements,1)),np.conjugate([inputs['dateN']]),axis=0)
+md.smb.pAir = np.append(np.tile(np.conjugate(inputs['pAir0']),(md.mesh.numberofelements,1)),np.conjugate([inputs['dateN']]),axis=0)
+md.smb.pAir = np.append(np.tile(np.conjugate(inputs['pAir0']),(md.mesh.numberofelements,1)),np.conjugate([inputs['dateN']]),axis=0)
+md.smb.Vz = np.tile(np.conjugate(inputs['LP']['Vz']),(md.mesh.numberofelements,1)).flatten()
+md.smb.Tz = np.tile(np.conjugate(inputs['LP']['Tz']),(md.mesh.numberofelements,1)).flatten()
+md.smb.Tmean = np.tile(np.conjugate(inputs['LP']['Tmean']),(md.mesh.numberofelements,1)).flatten()
+md.smb.C = np.tile(np.conjugate(inputs['LP']['C']),(md.mesh.numberofelements,1)).flatten()
+
+#smb settings
+md.smb.requested_outputs = ['SmbDz','SmbT','SmbD','SmbRe','SmbGdn','SmbGsp','SmbEC','SmbA','SmbMassBalance']
+
+#only run smb core: 
+md.transient.isstressbalance = 0
+md.transient.ismasstransport = 0
+md.transient.isthermal = 0
+
+#time stepping: 
+md.timestepping.start_time = 1965.
+md.timestepping.final_time = 1966.
+md.timestepping.time_step = 1. / 365.0
+md.timestepping.interp_forcings = 0.
+
+#Run transient
+md = solve(md,'Transient')
+
+#Fields and tolerances to track changes
+field_names      = ['SmbDz','SmbT' ,'SmbD' ,'SmbRe','SmbGdn','SmbGsp','SmbA' ,'SmbEC','SmbMassBalance']
+field_tolerances = [5e-4,5e-5,0.0006,0.0002,1e-5,0.0003,2e-5,2e-7,1e-7]
+#shape is different in python solution (fixed using reshape) which can cause test failure:
+field_values = [
+       md.results.TransientSolution[-1].SmbDz[0,0:240].reshape(1,-1),
+       md.results.TransientSolution[-1].SmbT[0,0:240].reshape(1,-1),
+       md.results.TransientSolution[-1].SmbD[0,0:240].reshape(1,-1),
+       md.results.TransientSolution[-1].SmbRe[0,0:240].reshape(1,-1),
+       md.results.TransientSolution[-1].SmbGdn[0,0:240].reshape(1,-1),
+       md.results.TransientSolution[-1].SmbGsp[0,0:240].reshape(1,-1),
+       md.results.TransientSolution[-1].SmbA[0,0:240].reshape(1,-1),
+       md.results.TransientSolution[-1].SmbEC[0],
+       md.results.TransientSolution[-1].SmbMassBalance[0],
+       ]
Index: ../trunk-jpl/src/m/solve/WriteData.py
===================================================================
--- ../trunk-jpl/src/m/solve/WriteData.py       (revision 22266)
+++ ../trunk-jpl/src/m/solve/WriteData.py       (revision 22267)
@@ -50,9 +50,9 @@
        if np.size(data) > 1 and np.size(data,0)==timeserieslength:
                yts = options.getfieldvalue('yts')
                if np.ndim(data) > 1:
-                       data[-1,:] = data[-1,:]*yts
+                       data[-1,:] = yts*data[-1,:]
                else:
-                       data[-1] = data[-1]*yts
+                       data[-1] = yts*data[-1]
 
        #Step 1: write the enum to identify this record uniquely
        fid.write(struct.pack('i',len(name)))
Index: ../trunk-jpl/src/m/solve/parseresultsfromdisk.py
===================================================================
--- ../trunk-jpl/src/m/solve/parseresultsfromdisk.py    (revision 22266)
+++ ../trunk-jpl/src/m/solve/parseresultsfromdisk.py    (revision 22267)
@@ -174,8 +174,6 @@
                yts=md.constants.yts
                if fieldname=='BalancethicknessThickeningRate':
                        field = field*yts
-               elif fieldname=='Time':
-                       field = field/yts
                elif fieldname=='HydrologyWaterVx':
                        field = field*yts
                elif fieldname=='HydrologyWaterVy':
@@ -190,17 +188,36 @@
                        field = field*yts
                elif fieldname=='BasalforcingsGroundediceMeltingRate':
                        field = field*yts
+               elif fieldname=='BasalforcingsFloatingiceMeltingRate':
+                       field = field*yts
                elif fieldname=='TotalFloatingBmb':
-                       field = field/10.**12.*yts #(GigaTon/year)
+                       field = field/10.**12*yts #(GigaTon/year)
                elif fieldname=='TotalGroundedBmb':
-                       field = field/10.**12.*yts #(GigaTon/year)
+                       field = field/10.**12*yts #(GigaTon/year)
                elif fieldname=='TotalSmb':
-                       field = field/10.**12.*yts #(GigaTon/year)
+                       field = field/10.**12*yts #(GigaTon/year)
                elif fieldname=='SmbMassBalance':
                        field = field*yts
+               elif fieldname=='SmbPrecipitation':
+                       field = field*yts
+               elif fieldname=='SmbRunoff':
+                       field = field*yts
+               elif fieldname=='SmbEC':
+                       field = field*yts
+               elif fieldname=='SmbAccumulation':
+                       field = field*yts
+               elif fieldname=='SmbMelt':
+                       field = field*yts
+               elif fieldname=='SmbDz_add':
+                       field = field*yts
+               elif fieldname=='SmbM_add':
+                       field = field*yts
                elif fieldname=='CalvingCalvingrate':
                        field = field*yts
 
+               if time !=-9999:
+                       time = time/yts
+
                saveres=OrderedDict()
                saveres['fieldname']=fieldname
                saveres['time']=time
Index: ../trunk-jpl/src/m/geometry/NowickiProfile.py
===================================================================
--- ../trunk-jpl/src/m/geometry/NowickiProfile.py       (nonexistent)
+++ ../trunk-jpl/src/m/geometry/NowickiProfile.py       (revision 22267)
@@ -0,0 +1,44 @@
+import numpy as np
+
+def NowickiProfile(x):
+       """
+       NOWICKIPROFILE - Create profile at the transition zone based on Sophie Nowicki's thesis
+
+       Usage:
+               [b h] = NowickiProfile(x)
+
+               - h = ice thickness
+               - b = ice base
+               - x = along flow coordinate
+       """
+       #Constant for theoretical profile
+       delta = 0.1             #ratio of water density and ice density -1
+       hg    = 1.              #ice thickness at grounding line
+       sea   = hg / (1+delta)  #sea level
+       lamda = 0.1             #ration of deviatoric stress and water pressure
+       beta  = 5.              #friction coefficient
+       ms    = 0.005           #surface accumulation rat
+       mu    = 5.              #viscosity
+       q     = 0.801           #ice mass flux
+
+       #mesh parameters
+       b = np.zeros((np.size(x),))
+       h = np.zeros((np.size(x),))
+       s = np.zeros((np.size(x),))
+
+       #upstream of the GL
+       for i in range(np.size(x) / 2):
+               ss = np.roots([1, 4 * lamda * beta, 0, 0, 6 * lamda * ms * x[i]**2 +
+                               12 * lamda * q * x[i] - hg**4 - 4 * lamda * beta * hg**3])
+               for j in range(4):
+                       if (np.isreal(ss[j]) > 0) and (np.imag(ss[j]) == 0):
+                               s[i] = ss[j]
+               h[i] = s[i]
+               b[i] = 0.
+
+       #downstream of the GL
+       for i in range(np.size(x) / 2, np.size(x)):
+               h[i] = (x[i] / (4. * (delta+1) * q) + hg**(-2))**(-0.5) # ice thickness for ice shelf from (3.1)
+               b[i] = sea - h[i] * (1. / (1+delta))
+
+       return [b, h, sea]
Index: ../trunk-jpl/src/m/consistency/checkfield.py
===================================================================
--- ../trunk-jpl/src/m/consistency/checkfield.py        (revision 22266)
+++ ../trunk-jpl/src/m/consistency/checkfield.py        (revision 22267)
@@ -83,7 +83,7 @@
        #Check numel
        if options.exist('numel'):
                fieldnumel=options.getfieldvalue('numel')
-               if np.size(field) not in fieldnumel:
+               if (type(fieldnumel) == int and np.size(field) != fieldnumel) or (type(fieldnumel) == list and np.size(field) not in fieldnumel):
                        if   len(fieldnumel)==1:
                                md = md.checkmessage(options.getfieldvalue('message',\
                                        "field '%s' size should be %d" % (fieldname,fieldnumel)))
@@ -100,6 +100,7 @@
                        md = md.checkmessage(options.getfieldvalue('message',\
                                "NaN values found in field '%s'" % fieldname))
 
+
        #check Inf
        if options.getfieldvalue('Inf',0):
                if np.any(np.isinf(field)):
@@ -106,6 +107,7 @@
                        md = md.checkmessage(options.getfieldvalue('message',\
                                "Inf values found in field '%s'" % fieldname))
 
+
        #check cell
        if options.getfieldvalue('cell',0):
                if not isinstance(field,(tuple,list,dict)):
@@ -128,13 +130,28 @@
 
        #check greater
        if options.exist('>='):
-               lowerbound=options.getfieldvalue('>=')
-               if np.any(field<lowerbound):
+               lowerbound = options.getfieldvalue('>=')
+               field2 = np.reshape(field,(np.prod(np.shape(field)),1),order='F').copy()
+               if options.getfieldvalue('timeseries',0):
+                       field2 = np.reshape(field[:-1],np.prod(np.shape(field[:-1])),1)
+
+               if options.getfieldvalue('singletimeseries',0):
+                       field2 = np.reshape(field[0],np.prod(np.shape(field[0])),1)
+
+               if np.any(field2<lowerbound):
                        md = md.checkmessage(options.getfieldvalue('message',\
                                "field '%s' should have values above %d" % (fieldname,lowerbound)))
+
        if options.exist('>'):
                lowerbound=options.getfieldvalue('>')
-               if np.any(field<=lowerbound):
+               field2 = np.reshape(field,(np.prod(np.shape(field)),1),order='F').copy()
+               if options.getfieldvalue('timeseries',0):
+                       field2 = np.reshape(field[:-1],np.prod(np.shape(field[:-1])),1)
+
+               if options.getfieldvalue('singletimeseries',0):
+                       field2 = np.reshape(field[0],np.prod(np.shape(field[0])),1)
+
+               if np.any(field2<=lowerbound):
                        md = md.checkmessage(options.getfieldvalue('message',\
                                "field '%s' should have values above %d" % (fieldname,lowerbound)))
 
@@ -141,12 +158,26 @@
        #check smaller
        if options.exist('<='):
                upperbound=options.getfieldvalue('<=')
-               if np.any(field>upperbound):
+               field2 = np.reshape(field,(np.prod(np.shape(field)),1),order='F').copy()
+               if options.getfieldvalue('timeseries',0):
+                       field2 = np.reshape(field[:-1],np.prod(np.shape(field[:-1])),1)
+
+               if options.getfieldvalue('singletimeseries',0):
+                       field2 = np.reshape(field[0],np.prod(np.shape(field[0])),1)
+
+               if np.any(field2>upperbound):
                        md = md.checkmessage(options.getfieldvalue('message',\
                                "field '%s' should have values below %d" % (fieldname,upperbound)))
        if options.exist('<'):
                upperbound=options.getfieldvalue('<')
-               if np.any(field>=upperbound):
+               field2 = np.reshape(field,(np.prod(np.shape(field)),1),order='F').copy()
+               if options.getfieldvalue('timeseries',0):
+                       field2 = np.reshape(field[:-1],np.prod(np.shape(field[:-1])),1)
+
+               if options.getfieldvalue('singletimeseries',0):
+                       field2 = np.reshape(field[0],np.prod(np.shape(field[0])),1)
+
+               if np.any(field2>=upperbound):
                        md = md.checkmessage(options.getfieldvalue('message',\
                                "field '%s' should have values below %d" % (fieldname,upperbound)))
 
@@ -163,15 +194,15 @@
 
        #Check forcings (size and times)
        if options.getfieldvalue('timeseries',0):
-               if   np.size(field,0)==md.mesh.numberofvertices:
+               if np.size(field,0)==md.mesh.numberofvertices or np.size(field,0)==md.mesh.numberofelements:
                        if np.ndim(field)>1 and not np.size(field,1)==1:
                                md = md.checkmessage(options.getfieldvalue('message',\
                                        "field '%s' should have only one column as there are md.mesh.numberofvertices lines" % fieldname))
-               elif np.size(field,0)==md.mesh.numberofvertices+1 or np.size(field,0)==2:
-                       if not all(field[-1,:]==np.sort(field[-1,:])):
+               elif np.size(field,0)==md.mesh.numberofvertices+1 or np.size(field,0)==md.mesh.numberofelements+1:
+                       if np.ndim(field) > 1 and not all(field[-1,:]==np.sort(field[-1,:])):
                                md = md.checkmessage(options.getfieldvalue('message',\
                                        "field '%s' columns should be sorted chronologically" % fieldname))
-                       if any(field[-1,0:-1]==field[-1,1:]):
+                       if np.ndim(field) > 1 and any(field[-1,0:-1]==field[-1,1:]):
                                md = md.checkmessage(options.getfieldvalue('message',\
                                        "field '%s' columns must not contain duplicate timesteps" % fieldname))
                else:
@@ -197,3 +228,4 @@
 
        return md
 
+
Index: ../trunk-jpl/src/m/mech/newforcing.py
===================================================================
--- ../trunk-jpl/src/m/mech/newforcing.py       (nonexistent)
+++ ../trunk-jpl/src/m/mech/newforcing.py       (revision 22267)
@@ -0,0 +1,34 @@
+import numpy as np
+
+def newforcing(t0,t1,deltaT,f0,f1,nodes):
+       '''
+FUNCTION NEWFORCING Build forcing that extends temporally from t0 to t1, and in magnitude from f0 to f1. Equal time 
+                    and magnitude spacing. 
+
+       Usage: forcing=newforcing(t0,t1,deltaT,f0,f1,nodes);  
+       Where: 
+          t0:t1: time interval. 
+          deltaT: time step
+          f0:f1: magnitude interval.
+          nodes: number of vertices where we have a temporal forcing
+
+       Example: 
+           md.smb.mass_balance=newforcing(md.timestepping.start_time,md.timestepping.final_time,
+                                          md.timestepping.time_step,-1,+2,md.mesh.numberofvertices)
+       '''
+       #Number of time steps: 
+       nsteps = (t1 - t0) / deltaT + 1
+
+       #delta forcing:
+       deltaf = (f1 - f0) / (nsteps - 1)
+
+       #creates times:
+       times = np.arange(t0,t1+deltaT,deltaT)  #Add deltaT to fix python/matlab discrepency
+
+       #create forcing:
+       forcing = np.arange(f0,f1+deltaf,deltaf)#Add deltaf to fix python/matlab discrepency
+
+       #replicate for all nodes
+       forcing = np.tile(forcing, (nodes+1,1))
+       forcing[-1,:] = times
+       return forcing
Index: ../trunk-jpl/src/m/dev/devpath.py
===================================================================
--- ../trunk-jpl/src/m/dev/devpath.py   (revision 22266)
+++ ../trunk-jpl/src/m/dev/devpath.py   (revision 22267)
@@ -30,6 +30,7 @@
 
 #Manual imports for commonly used functions
 from plotmodel import plotmodel
+from runme import runme
 
 #c = get_ipython().config
 #c.InteractiveShellApp.exec_lines = []
Index: ../trunk-jpl/src/m/print/printmodel.py
===================================================================
--- ../trunk-jpl/src/m/print/printmodel.py      (nonexistent)
+++ ../trunk-jpl/src/m/print/printmodel.py      (revision 22267)
@@ -0,0 +1,109 @@
+
+import numpy as np
+#from model import *
+#from socket import gethostname
+#from bamg import *
+#from setmask import *
+#from parameterize import *
+#from setflowequation import *
+#from solve import *
+from pairoptions import *
+
+def printmodel(filename,format,*args):
+'''    PRINTMODEL - save an image of current figure
+
+       filename: output name of image file (no extension)
+       format: image format (ex: 'tiff','jpg','pdf') 
+
+       List of options to printfmodel: 
+
+       figure: number of figure to print (default: current figure)
+       resolution: use higher resolution to anti-alias (default 2)
+       margin: add margin around final image  
+       marginsize: size of margin around final image (default 5)
+       frame: add frame around final image
+       framesize: size of frame around final image (default 5)
+       framecolor: color of frame around final image (default 'black')
+       trim: trim empty space around image (default 'off')
+       hardcopy: 'off' to impose MATLAB to use the same colors (default 'off')
+   
+       Usage:
+               printmodel(filename,format,varargin);
+
+       Examples:
+               printmodel('image','tiff')
+               printmodel('image','eps','margin','on','frame','on','hardcopy','on')
+'''
+
+#get options: 
+options = pairoptions(*args)
+
+#set defaults
+options = addfielddefault(options,'figure','gcf')
+options = addfielddefault(options,'format','tiff')
+options = addfielddefault(options,'resolution',1)
+options = addfielddefault(options,'margin','on')
+options = addfielddefault(options,'marginsize',25)
+options = addfielddefault(options,'frame','on')
+options = addfielddefault(options,'framesize',3)
+options = addfielddefault(options,'framecolor','black')
+options = addfielddefault(options,'trim','on')
+options = addfielddefault(options,'hardcopy','off')
+
+#get fig: 
+fig = getfieldvalue(options,'figure')
+if len(fig) == 1:
+       fig=gcf
+else:
+       figure(fig)
+       fig=gcf
+
+#In auto mode, MATLAB prints the figure the same size as it appears on the computer screen, centered on the page
+set(fig, 'PaperPositionMode', 'auto');
+
+#InvertHardcopy off imposes MATLAB to use the same colors
+set(fig, 'InvertHardcopy', getfieldvalue(options,'hardcopy'));
+
+#we could have several formats, as a cell array of strings.
+formats=format;
+if ~iscell(formats),
+       formats={formats};
+end
+
+#loop on formats:
+for i=1:length(formats),
+       format=formats{i};
+
+       #Use higher resolution to anti-alias and use zbuffer to have smooth colors
+       print(fig, '-zbuffer','-dtiff',['-r' num2str(get(0,'ScreenPixelsPerInch')*getfieldvalue(options,'resolution'))],filename);
+
+       #some trimming involved? 
+       if ~strcmpi(format,'pdf'),
+               if strcmpi(getfieldvalue(options,'trim'),'on'),
+                       system(['convert -trim ' filename '.tif ' filename '.tif']);
+               end
+       end
+
+       #margin?
+       if ~strcmpi(format,'pdf'),
+               if strcmpi(getfieldvalue(options,'margin'),'on'),
+                       marginsize=getfieldvalue(options,'marginsize');
+                       system(['convert -border ' num2str(marginsize) 'x' num2str(marginsize) ' -bordercolor "white" ' filename '.tif ' filename '.tif']);
+               end
+       end
+
+       #frame?
+       if ~strcmpi(format,'pdf'),
+               if strcmpi(getfieldvalue(options,'frame'),'on'),
+                       framesize=getfieldvalue(options,'framesize');
+                       framecolor=getfieldvalue(options,'framecolor');
+                       system(['convert -border ' num2str(framesize) 'x' num2str(framesize) ' -bordercolor "' framecolor '" ' filename '.tif ' filename '.tif']);
+               end
+       end
+
+       #convert image to correct format
+       if ~strcmpi(format,'tiff') & ~strcmpi(format,'tif'),
+               system(['convert ' filename '.tif ' filename '.' format]);
+               delete([ filename '.tif']);
+       end
+end
Index: ../trunk-jpl/src/m/classes/thermal.py
===================================================================
--- ../trunk-jpl/src/m/classes/thermal.py       (revision 22266)
+++ ../trunk-jpl/src/m/classes/thermal.py       (revision 22267)
@@ -100,14 +100,18 @@
                md = checkfield(md,'fieldname','thermal.requested_outputs','stringrow',1)
 
                if 'EnthalpyAnalysis' in analyses and md.thermal.isenthalpy and md.mesh.dimension()==3:
-                       pos=np.where(~np.isnan(md.thermal.spctemperature[0:md.mesh.numberofvertices]))
+                       TEMP = md.thermal.spctemperature[:-1].flatten(-1)
+                       pos=np.where(~np.isnan(TEMP))
                        try:
                                spccol=np.size(md.thermal.spctemperature,1)
                        except IndexError:
                                spccol=1
-                       replicate=np.tile(md.geometry.surface-md.mesh.z,(spccol))
-                       control=md.materials.meltingpoint-md.materials.beta*md.materials.rho_ice*md.constants.g*replicate
-                       md = checkfield(md,'fieldname','thermal.spctemperature','field',md.thermal.spctemperature[pos],'<=',control[pos],'message',"spctemperature should be below the adjusted melting point")
+
+                       replicate=np.tile(md.geometry.surface-md.mesh.z,(spccol)).flatten(-1)
+
+                       control=md.materials.meltingpoint-md.materials.beta*md.materials.rho_ice*md.constants.g*replicate+10**-5
+
+                       md = checkfield(md,'fieldname','thermal.spctemperature','field',md.thermal.spctemperature.flatten(-1)[pos],'<=',control[pos],'message',"spctemperature should be below the adjusted melting point")
                        md = checkfield(md,'fieldname','thermal.isenthalpy','numel',[1],'values',[0,1])
                        md = checkfield(md,'fieldname','thermal.isdynamicbasalspc','numel',[1],'values',[0,1]);
                        if(md.thermal.isenthalpy):
Index: ../trunk-jpl/src/m/classes/calvingminthickness.py
===================================================================
--- ../trunk-jpl/src/m/classes/calvingminthickness.py   (nonexistent)
+++ ../trunk-jpl/src/m/classes/calvingminthickness.py   (revision 22267)
@@ -0,0 +1,52 @@
+from fielddisplay import fielddisplay
+from checkfield import checkfield
+from WriteData import WriteData
+
+class calvingminthickness(object):
+       """
+       CALVINGMINTHICKNESS class definition
+
+          Usage:
+             calvingminthickness=calvingminthickness()
+       """
+
+       def __init__(self): # {{{
+
+               self.min_thickness = 0.
+               self.meltingrate   = float('NaN')
+
+               #set defaults
+               self.setdefaultparameters()
+
+       #}}}
+       def __repr__(self): # {{{
+               string='   Calving Minimum thickness:'
+               string="%s\n%s"%(string,fielddisplay(self,'min_thickness','minimum thickness below which no ice is allowed'))
+               string="%s\n%s"%(string,fielddisplay(self,'meltingrate','melting rate at given location [m/a]'))
+               return string
+       #}}}
+       def extrude(self,md): # {{{
+               self.meltingrate=project3d(md,'vector',self.meltingrate,'type','node')
+               return self
+       #}}}
+       def setdefaultparameters(self): # {{{
+
+               #minimum thickness is 100 m by default
+               self.min_thickness = 100.
+       #}}}
+       def checkconsistency(self,md,solution,analyses):    # {{{
+
+               #Early return
+               if solution == 'TransientSolution' or md.transient.ismovingfront == 0:
+                       return
+
+               md = checkfield(md,'fieldname','calving.min_thickness','>',0,'NaN',1,'Inf',1)
+               md = checkfield(md,'fieldname','calving.meltingrate','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices],'>=',0)
+               return md
+       # }}}
+       def marshall(self,prefix,md,fid):    # {{{
+               yts=md.constants.yts
+               WriteData(fid,prefix,'name','md.calving.law','data',4,'format','Integer')
+               WriteData(fid,prefix,'object',self,'fieldname','min_thickness','format','Double')
+               WriteData(fid,prefix,'object',self,'fieldname','meltingrate','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'scale',1./yts)
+       # }}}
Index: ../trunk-jpl/src/m/classes/calvingdev.py
===================================================================
--- ../trunk-jpl/src/m/classes/calvingdev.py    (nonexistent)
+++ ../trunk-jpl/src/m/classes/calvingdev.py    (revision 22267)
@@ -0,0 +1,60 @@
+from fielddisplay import fielddisplay
+from project3d import project3d
+from checkfield import checkfield
+from WriteData import WriteData
+
+class calvingdev(object):
+       """
+       CALVINGDEV class definition
+
+          Usage:
+             calvingdev=calvingdev();
+       """
+
+       def __init__(self): # {{{
+
+               self.stress_threshold_groundedice = 0.
+               self.stress_threshold_floatingice = 0.
+               self.meltingrate   = float('NaN')
+
+               #set defaults
+               self.setdefaultparameters()
+
+       #}}}
+       def __repr__(self): # {{{
+               string='   Calving Pi parameters:'
+               string="%s\n%s"%(string,fielddisplay(self,'stress_threshold_groundedice','sigma_max applied to grounded ice only [Pa]'))
+               string="%s\n%s"%(string,fielddisplay(self,'stress_threshold_floatingice','sigma_max applied to floating ice only [Pa]'))
+
+               string="%s\n%s"%(string,fielddisplay(self,'meltingrate','melting rate at given location [m/a]'))
+               return string
+       #}}}
+       def extrude(self,md): # {{{
+               self.meltingrate=project3d(md,'vector',self.meltingrate,'type','node')
+               return self
+       #}}}
+       def setdefaultparameters(self): # {{{
+               #Default sigma max
+               self.stress_threshold_groundedice = 1e6
+               self.stress_threshold_floatingice = 150e3
+               return self
+       #}}}
+       def checkconsistency(self,md,solution,analyses):    # {{{
+               #Early return
+               if solution == 'TransientSolution' or md.transient.ismovingfront == 0:
+                       return
+
+               md = checkfield(md,'fieldname','calving.stress_threshold_groundedice','>',0,'nan',1,'Inf',1)
+               md = checkfield(md,'fieldname','calving.stress_threshold_floatingice','>',0,'nan',1,'Inf',1)
+               md = checkfield(md,'fieldname','calving.meltingrate','NaN',1,'Inf',1,'timeseries',1,'>=',0)
+
+               return md
+       # }}}
+       def marshall(self,prefix,md,fid):    # {{{
+               yts=md.constants.yts
+
+               WriteData(fid,prefix,'name','md.calving.law','data',2,'format','Integer')
+               WriteData(fid,prefix,'object',self,'fieldname','stress_threshold_groundedice','format','DoubleMat','mattype',1)
+               WriteData(fid,prefix,'object',self,'fieldname','stress_threshold_floatingice','format','DoubleMat','mattype',1)
+               WriteData(fid,prefix,'object',self,'fieldname','meltingrate','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts,'scale',1./yts)
+       # }}}
Index: ../trunk-jpl/src/m/classes/mismipbasalforcings.py
===================================================================
--- ../trunk-jpl/src/m/classes/mismipbasalforcings.py   (revision 22266)
+++ ../trunk-jpl/src/m/classes/mismipbasalforcings.py   (revision 22267)
@@ -39,8 +39,11 @@
     #}}}
     def initialize(self,md): # {{{
         if np.all(np.isnan(self.groundedice_melting_rate)):
-            self.groundedice_melting_rate=np.zeros(md.mesh.numberofvertices)
+            self.groundedice_melting_rate=np.zeros((md.mesh.numberofvertices))
             print ' no basalforcings.groundedice_melting_rate specified: values set as zero'
+       if np.all(np.isnan(self.geothermalflux)):
+                       self.geothermalflux=np.zeros((md.mesh.numberofvertices))
+                       print "      no basalforcings.geothermalflux specified: values set as zero"
         return self
     #}}}
     def setdefaultparameters(self): # {{{
@@ -83,7 +86,7 @@
             print 'WARNING: value of yts for MISMIP+ runs different from ISSM default!'
 
         floatingice_melting_rate = np.zeros((md.mesh.numberofvertices))
-        floatingice_melting_rate = md.basalforcings.meltrate_factor*np.tanh((md.geometry.base-md.geometry.bed)/md.basalforcings.threshold_thickness)*np.amax(md.basalforcings.upperdepth_melt-md.geometry.base,0)
+        floatingice_melting_rate = md.basalforcings.meltrate_factor*np.tanh((md.geometry.base-md.geometry.bed)/md.basalforcings.threshold_thickness)*(md.basalforcings.upperdepth_melt-md.geometry.base)
 
        WriteData(fid,prefix,'name','md.basalforcings.model','data',3,'format','Integer')
        WriteData(fid,prefix,'data',floatingice_melting_rate,'format','DoubleMat','name','md.basalforcings.floatingice_melting_rate','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts)
Index: ../trunk-jpl/src/m/classes/SMBgradientsela.py
===================================================================
--- ../trunk-jpl/src/m/classes/SMBgradientsela.py       (revision 22266)
+++ ../trunk-jpl/src/m/classes/SMBgradientsela.py       (revision 22267)
@@ -8,7 +8,7 @@
        SMBgradientsela Class definition
 
           Usage:
-             SMBgradientsela=SMBgradientsela();
+             SMBgradientsela=SMBgradientsela()
        """
 
        def __init__(self): # {{{
@@ -18,11 +18,12 @@
                self.b_max   = float('NaN')
                self.b_min   = float('NaN')
                self.requested_outputs      = []
+               self.setdefaultparameters()
                #}}}
        def __repr__(self): # {{{
-               string="   surface forcings parameters:"
+               string = "   surface forcings parameters:"
+               string+= '\n   SMB gradients ela parameters:'
 
-               string="%s\n%s"%(string,fielddisplay(self,'issmbgradientsela','is smb gradients ela method activated (0 or 1, default is 0)'))
                string="%s\n%s"%(string,fielddisplay(self,'ela',' equilibrium line altitude from which deviation is used to calculate smb using the smb gradients ela method [m a.s.l.]'))
                string="%s\n%s"%(string,fielddisplay(self,'b_pos',' vertical smb gradient (dB/dz) above ela'))
                string="%s\n%s"%(string,fielddisplay(self,'b_neg',' vertical smb gradient (dB/dz) below ela'))
@@ -46,6 +47,11 @@
 
                return self
        #}}}
+       def setdefaultparameters(self): # {{{
+               self.b_max=9999.
+               self.b_min=-9999.
+               return self
+       #}}}
        def checkconsistency(self,md,solution,analyses):    # {{{
 
                if 'MasstransportAnalysis' in analyses:
@@ -55,7 +61,7 @@
                        md = checkfield(md,'fieldname','smb.b_max','timeseries',1,'NaN',1,'Inf',1)
                        md = checkfield(md,'fieldname','smb.b_min','timeseries',1,'NaN',1,'Inf',1)
 
-               md = checkfield(md,'fieldname','masstransport.requested_outputs','stringrow',1)
+               md = checkfield(md,'fieldname','smb.requested_outputs','stringrow',1)
                return md
        # }}}
        def marshall(self,prefix,md,fid):    # {{{
Index: ../trunk-jpl/src/m/classes/constants.py
===================================================================
--- ../trunk-jpl/src/m/classes/constants.py     (revision 22266)
+++ ../trunk-jpl/src/m/classes/constants.py     (revision 22267)
@@ -11,10 +11,10 @@
        """
 
        def __init__(self): # {{{
-               self.g                    = 0
-               self.omega                = 0
-               self.yts                  = 0
-               self.referencetemperature = 0
+               self.g                    = 0.
+               self.omega                = 0.
+               self.yts                  = 0.
+               self.referencetemperature = 0.
                
                #set defaults
                self.setdefaultparameters()
@@ -48,10 +48,10 @@
        #}}}
        def checkconsistency(self,md,solution,analyses):    # {{{
 
-               md = checkfield(md,'fieldname','constants.g','>',0,'size',[1])
-               md = checkfield(md,'fieldname','constants.omega','>=',0,'size',[1])
-               md = checkfield(md,'fieldname','constants.yts','>',0,'size',[1])
-               md = checkfield(md,'fieldname','constants.referencetemperature','size',[1])
+               md = checkfield(md,'fieldname','constants.g','>=',0,'size',[1,1])
+               md = checkfield(md,'fieldname','constants.omega','>=',0,'size',[1,1])
+               md = checkfield(md,'fieldname','constants.yts','>',0,'size',[1,1])
+               md = checkfield(md,'fieldname','constants.referencetemperature','size',[1,1])
 
                return md
        # }}}
Index: ../trunk-jpl/src/m/classes/hydrologysommers.py
===================================================================
--- ../trunk-jpl/src/m/classes/hydrologysommers.py      (nonexistent)
+++ ../trunk-jpl/src/m/classes/hydrologysommers.py      (revision 22267)
@@ -0,0 +1,90 @@
+from fielddisplay import fielddisplay
+from checkfield import checkfield
+from WriteData import WriteData
+
+class hydrologysommers(object):
+       """
+       HYDROLOGYSOMMERS class definition
+
+          Usage:
+             hydrologysommers=hydrologysommers()
+       """
+
+       def __init__(self): # {{{
+               self.head            = float('NaN')
+               self.gap_height      = float('NaN')
+               self.bump_spacing    = float('NaN')
+               self.bump_height     = float('NaN')
+               self.englacial_input = float('NaN')
+               self.moulin_input    = float('NaN')
+               self.reynolds        = float('NaN')
+               self.spchead         = float('NaN')
+               self.neumannflux     = float('NaN')
+               self.relaxation      = 0
+               self.storage         = 0
+
+               #set defaults
+               self.setdefaultparameters()
+
+               #}}}
+       def __repr__(self): # {{{
+               
+               string='   hydrologysommers solution parameters:'
+               string="%s\n%s"%(string,fielddisplay(self,'head','subglacial hydrology water head (m)'))
+               string="%s\n%s"%(string,fielddisplay(self,'gap_height','height of gap separating ice to bed (m)'))
+               string="%s\n%s"%(string,fielddisplay(self,'bump_spacing','characteristic bedrock bump spacing (m)'))
+               string="%s\n%s"%(string,fielddisplay(self,'bump_height','characteristic bedrock bump height (m)'))
+               string="%s\n%s"%(string,fielddisplay(self,'englacial_input','liquid water input from englacial to subglacial system (m/yr)'))
+               string="%s\n%s"%(string,fielddisplay(self,'moulin_input','liquid water input from moulins (at the vertices) to subglacial system (m^3/s)'))
+               string="%s\n%s"%(string,fielddisplay(self,'reynolds','Reynolds'' number'))
+               string="%s\n%s"%(string,fielddisplay(self,'neumannflux','water flux applied along the model boundary (m^2/s)'))
+               string="%s\n%s"%(string,fielddisplay(self,'spchead','water head constraints (NaN means no constraint) (m)'))
+               string="%s\n%s"%(string,fielddisplay(self,'relaxation','under-relaxation coefficient for nonlinear iteration'))
+               string="%s\n%s"%(string,fielddisplay(self,'storage','englacial storage coefficient (void ratio)'))
+               return string
+               #}}}
+       def extrude(self,md): # {{{
+               return self
+       #}}}
+       def setdefaultparameters(self): # {{{
+               # Set under-relaxation parameter to be 1 (no under-relaxation of nonlinear iteration)   
+               self.relaxation=1
+               self.storage=0
+               return self
+       #}}}
+       def checkconsistency(self,md,solution,analyses):    # {{{
+               
+               #Early return
+               if 'HydrologySommersAnalysis' not in analyses:
+                       return md
+
+               md = checkfield(md,'fieldname','hydrology.head','size',[md.mesh.numberofvertices],'NaN',1,'Inf',1)
+               md = checkfield(md,'fieldname','hydrology.gap_height','>=',0,'size',[md.mesh.numberofelements],'NaN',1,'Inf',1)
+               md = checkfield(md,'fieldname','hydrology.bump_spacing','>',0,'size',[md.mesh.numberofelements],'NaN',1,'Inf',1)
+               md = checkfield(md,'fieldname','hydrology.bump_height','>=',0,'size',[md.mesh.numberofelements],'NaN',1,'Inf',1)
+               md = checkfield(md,'fieldname','hydrology.englacial_input','>=',0,'NaN',1,'Inf',1,'timeseries',1)
+               md = checkfield(md,'fieldname','hydrology.moulin_input','>=',0,'NaN',1,'Inf',1,'timeseries',1)
+               md = checkfield(md,'fieldname','hydrology.reynolds','>',0,'size',[md.mesh.numberofelements],'NaN',1,'Inf',1)
+               md = checkfield(md,'fieldname','hydrology.neumannflux','timeseries',1,'NaN',1,'Inf',1)
+               md = checkfield(md,'fieldname','hydrology.spchead','size',[md.mesh.numberofvertices])
+               md = checkfield(md,'fieldname','hydrology.relaxation','>=',0)   
+               md = checkfield(md,'fieldname','hydrology.storage','>=',0)
+
+               return md
+       # }}}
+       def marshall(self,prefix,md,fid):    # {{{
+               yts=md.constants.yts
+
+               WriteData(fid,prefix,'name','md.hydrology.model','data',3,'format','Integer')
+               WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','head','format','DoubleMat','mattype',1)
+               WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','gap_height','format','DoubleMat','mattype',2)
+               WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','bump_spacing','format','DoubleMat','mattype',2)
+               WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','bump_height','format','DoubleMat','mattype',2)
+               WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','englacial_input','format','DoubleMat','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts)
+               WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','moulin_input','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts)
+               WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','reynolds','format','DoubleMat','mattype',2)
+               WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','neumannflux','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts)
+               WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','spchead','format','DoubleMat','mattype',1)
+               WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','relaxation','format','Double')
+               WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','storage','format','Double')
+       # }}}
Index: ../trunk-jpl/src/m/classes/flowequation.py
===================================================================
--- ../trunk-jpl/src/m/classes/flowequation.py  (revision 22266)
+++ ../trunk-jpl/src/m/classes/flowequation.py  (revision 22267)
@@ -91,7 +91,7 @@
                md = checkfield(md,'fieldname','flowequation.isFS','numel',[1],'values',[0,1])
                md = checkfield(md,'fieldname','flowequation.fe_SSA','values',['P1','P1bubble','P1bubblecondensed','P2','P2bubble'])
                md = checkfield(md,'fieldname','flowequation.fe_HO' ,'values',['P1','P1bubble','P1bubblecondensed','P1xP2','P2xP1','P2','P2bubble','P1xP3','P2xP4'])
-               md = checkfield(md,'fieldname','flowequation.fe_FS' ,'values',['P1P1','P1P1GLS','MINIcondensed','MINI','TaylorHood','XTaylorHood','OneLayerP4z','CrouzeixRaviart'])
+               md = checkfield(md,'fieldname','flowequation.fe_FS' ,'values',['P1P1','P1P1GLS','MINIcondensed','MINI','TaylorHood','LATaylorHood','XTaylorHood','OneLayerP4z','CrouzeixRaviart','LACrouzeixRaviart'])
                md = checkfield(md,'fieldname','flowequation.borderSSA','size',[md.mesh.numberofvertices],'values',[0,1])
                md = checkfield(md,'fieldname','flowequation.borderHO','size',[md.mesh.numberofvertices],'values',[0,1])
                md = checkfield(md,'fieldname','flowequation.borderFS','size',[md.mesh.numberofvertices],'values',[0,1])
@@ -103,6 +103,9 @@
                if m.strcmp(md.mesh.domaintype(),'2Dhorizontal'):
                        md = checkfield(md,'fieldname','flowequation.vertex_equation','size',[md.mesh.numberofvertices],'values',[1,2])
                        md = checkfield(md,'fieldname','flowequation.element_equation','size',[md.mesh.numberofelements],'values',[1,2])
+               elif m.strcmp(md.mesh.domaintype(),'2Dvertical'):
+                       md = checkfield(md,'fieldname','flowequation.vertex_equation','size',[md.mesh.numberofvertices],'values',[2,4,5])
+                       md = checkfield(md,'fieldname','flowequation.element_equation','size',[md.mesh.numberofelements],'values',[2,4,5])
                elif m.strcmp(md.mesh.domaintype(),'3D'):
                        md = checkfield(md,'fieldname','flowequation.vertex_equation','size',[md.mesh.numberofvertices],'values',np.arange(0,8+1))
                        md = checkfield(md,'fieldname','flowequation.element_equation','size',[md.mesh.numberofelements],'values',np.arange(0,8+1))
Index: ../trunk-jpl/src/m/classes/basalforcings.py
===================================================================
--- ../trunk-jpl/src/m/classes/basalforcings.py (revision 22266)
+++ ../trunk-jpl/src/m/classes/basalforcings.py (revision 22267)
@@ -44,6 +44,9 @@
                if np.all(np.isnan(self.floatingice_melting_rate)):
                        self.floatingice_melting_rate=np.zeros((md.mesh.numberofvertices))
                        print "      no basalforcings.floatingice_melting_rate specified: values set as zero"
+               #if np.all(np.isnan(self.geothermalflux)):
+                       #self.geothermalflux=np.zeros((md.mesh.numberofvertices))
+                       #print "      no basalforcings.geothermalflux specified: values set as zero"
 
                return self
        #}}}
Index: ../trunk-jpl/src/m/classes/calvingvonmises.py
===================================================================
--- ../trunk-jpl/src/m/classes/calvingvonmises.py       (nonexistent)
+++ ../trunk-jpl/src/m/classes/calvingvonmises.py       (revision 22267)
@@ -0,0 +1,60 @@
+from fielddisplay import fielddisplay
+from project3d import project3d
+from checkfield import checkfield
+from WriteData import WriteData
+
+class calvingvonmises(object):
+       """
+       CALVINGVONMISES class definition
+
+          Usage:
+             calvingvonmises=calvingvonmises()
+       """
+
+       def __init__(self): # {{{
+
+               self.stress_threshold_groundedice = 0.
+               self.stress_threshold_floatingice = 0.
+               self.meltingrate   = float('NaN')
+
+               #set defaults
+               self.setdefaultparameters()
+
+       #}}}
+       def __repr__(self): # {{{
+               string='   Calving VonMises parameters:'
+               string="%s\n%s"%(string,fielddisplay(self,'stress_threshold_groundedice','sigma_max applied to grounded ice only [Pa]'))
+               string="%s\n%s"%(string,fielddisplay(self,'stress_threshold_floatingice','sigma_max applied to floating ice only [Pa]'))
+
+               string="%s\n%s"%(string,fielddisplay(self,'meltingrate','melting rate at given location [m/a]'))
+               return string
+       #}}}
+       def extrude(self,md): # {{{
+               self.meltingrate=project3d(md,'vector',self.meltingrate,'type','node')
+               return self
+       #}}}
+       def setdefaultparameters(self): # {{{
+               #Default sigma max
+               self.stress_threshold_groundedice = 1e6
+               self.stress_threshold_floatingice = 150e3
+               return self
+       #}}}
+       def checkconsistency(self,md,solution,analyses):    # {{{
+               #Early return
+               if solution == 'TransientSolution' or md.transient.ismovingfront == 0:
+                       return
+
+               md = checkfield(md,'fieldname','calving.stress_threshold_groundedice','>',0,'nan',1,'Inf',1)
+               md = checkfield(md,'fieldname','calving.stress_threshold_floatingice','>',0,'nan',1,'Inf',1)
+               md = checkfield(md,'fieldname','calving.meltingrate','NaN',1,'Inf',1,'timeseries',1,'>=',0)
+
+               return md
+       # }}}
+       def marshall(self,prefix,md,fid):    # {{{
+               yts=md.constants.yts
+
+               WriteData(fid,prefix,'name','md.calving.law','data',2,'format','Integer')
+               WriteData(fid,prefix,'object',self,'fieldname','stress_threshold_groundedice','format','DoubleMat','mattype',1)
+               WriteData(fid,prefix,'object',self,'fieldname','stress_threshold_floatingice','format','DoubleMat','mattype',1)
+               WriteData(fid,prefix,'object',self,'fieldname','meltingrate','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts,'scale',1./yts)
+       # }}}
Index: ../trunk-jpl/src/m/classes/transient.py
===================================================================
--- ../trunk-jpl/src/m/classes/transient.py     (revision 22266)
+++ ../trunk-jpl/src/m/classes/transient.py     (revision 22267)
@@ -11,7 +11,7 @@
        """
 
        def __init__(self): # {{{
-               self.issmb   = False
+               self.issmb             = False
                self.ismasstransport   = False
                self.isstressbalance   = False
                self.isthermal         = False
@@ -22,9 +22,9 @@
                self.ismovingfront     = False
                self.ishydrology       = False
                self.isslr             = False
+               self.iscoupler         = False
+               self.amr_frequency     = 0
                self.isoceancoupling   = False
-               self.iscoupler         = False
-               amr_frequency                     = 0
                self.requested_outputs = []
 
                #set defaults
@@ -80,6 +80,26 @@
                self.requested_outputs=[]
                return self
        #}}}
+       def deactivateall(self):#{{{
+               self.issmb             = False
+               self.ismasstransport   = False
+               self.isstressbalance   = False
+               self.isthermal         = False
+               self.isgroundingline   = False
+               self.isgia             = False
+               self.isesa             = False
+               self.isdamageevolution = False
+               self.ismovingfront     = False
+               self.ishydrology       = False
+               self.isslr             = False
+               self.isoceancoupling   = False
+               self.iscoupler         = False
+               self.amr_frequency     = 0
+
+               #default output
+               self.requested_outputs=[]
+               return self
+       #}}}
        def setdefaultparameters(self): # {{{
                
                #full analysis: Stressbalance, Masstransport and Thermal but no groundingline migration for now
Index: ../trunk-jpl/src/m/classes/mesh2dvertical.py
===================================================================
--- ../trunk-jpl/src/m/classes/mesh2dvertical.py        (nonexistent)
+++ ../trunk-jpl/src/m/classes/mesh2dvertical.py        (revision 22267)
@@ -0,0 +1,130 @@
+import numpy as np
+from fielddisplay import fielddisplay
+from checkfield import checkfield
+import MatlabFuncs as m
+from WriteData import WriteData
+
+class mesh2dvertical(object):
+       """
+       MESH2DVERTICAL class definition
+
+          Usage:
+             mesh2dvertical=mesh2dvertical();
+       """
+
+       def __init__(self): # {{{
+               self.x                           = float('NaN')
+               self.y                           = float('NaN')
+               self.elements                    = float('NaN')
+               self.numberofelements            = 0
+               self.numberofvertices            = 0
+               self.numberofedges               = 0
+               
+               self.lat                         = float('NaN')
+               self.long                        = float('NaN')
+               self.epsg                        = float('NaN')
+
+               self.vertexonboundary            = float('NaN')
+               self.vertexonbase                = float('NaN')
+               self.vertexonsurface             = float('NaN')
+
+               self.edges                       = float('NaN')
+               self.segments                    = float('NaN')
+               self.segmentmarkers              = float('NaN')
+               self.vertexconnectivity          = float('NaN')
+               self.elementconnectivity         = float('NaN')
+               self.average_vertex_connectivity = 0
+
+               #set defaults
+               self.setdefaultparameters()
+
+               #}}}
+       def __repr__(self): # {{{
+               string="   2D tria Mesh (vertical):" 
+
+               string="%s\n%s"%(string,"\n      Elements and vertices:")
+               string="%s\n%s"%(string,fielddisplay(self,"numberofelements","number of elements"))
+               string="%s\n%s"%(string,fielddisplay(self,"numberofvertices","number of vertices"))
+               string="%s\n%s"%(string,fielddisplay(self,"elements","vertex indices of the mesh elements"))
+               string="%s\n%s"%(string,fielddisplay(self,"x","vertices x coordinate [m]"))
+               string="%s\n%s"%(string,fielddisplay(self,"y","vertices y coordinate [m]"))
+               string="%s\n%s"%(string,fielddisplay(self,"edges","edges of the 2d mesh (vertex1 vertex2 element1 element2)"))
+               string="%s\n%s"%(string,fielddisplay(self,"numberofedges","number of edges of the 2d mesh"))
+
+               string="%s%s"%(string,"\n\n      Properties:")
+               string="%s\n%s"%(string,fielddisplay(self,"vertexonboundary","vertices on the boundary of the domain flag list"))
+               string="%s\n%s"%(string,fielddisplay(self,'vertexonbase','vertices on the bed of the domain flag list'))
+               string="%s\n%s"%(string,fielddisplay(self,'vertexonsurface','vertices on the surface of the domain flag list'))
+               string="%s\n%s"%(string,fielddisplay(self,"segments","edges on domain boundary (vertex1 vertex2 element)"))
+               string="%s\n%s"%(string,fielddisplay(self,"segmentmarkers","number associated to each segment"))
+               string="%s\n%s"%(string,fielddisplay(self,"vertexconnectivity","list of vertices connected to vertex_i"))
+               string="%s\n%s"%(string,fielddisplay(self,"elementconnectivity","list of vertices connected to element_i"))
+               string="%s\n%s"%(string,fielddisplay(self,"average_vertex_connectivity","average number of vertices connected to one vertex"))
+
+               string="%s%s"%(string,"\n\n      Projection:")
+               string="%s\n%s"%(string,fielddisplay(self,"lat","vertices latitude [degrees]"))
+               string="%s\n%s"%(string,fielddisplay(self,"long","vertices longitude [degrees]"))
+               string="%s\n%s"%(string,fielddisplay(self,"epsg","EPSG code (ex: 3413 for UPS Greenland, 3031 for UPS Antarctica)"))
+               return string
+               #}}}
+       def setdefaultparameters(self): # {{{
+               
+               #the connectivity is the averaged number of nodes linked to a
+               #given node through an edge. This connectivity is used to initially
+               #allocate memory to the stiffness matrix. A value of 16 seems to
+               #give a good memory/time ration. This value can be checked in
+               #trunk/test/Miscellaneous/runme.m
+               self.average_vertex_connectivity=25.
+
+               return self
+       #}}}
+       def checkconsistency(self,md,solution,analyses):    # {{{
+               if(solution=='LoveSolution'):
+                       return
+
+               md = checkfield(md,'fieldname','mesh.x','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices])
+               md = checkfield(md,'fieldname','mesh.y','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices])
+               md = checkfield(md,'fieldname','mesh.elements','NaN',1,'Inf',1,'>',0,'values',np.arange(1,md.mesh.numberofvertices+1))
+               md = checkfield(md,'fieldname','mesh.elements','size',[md.mesh.numberofelements,3])
+               if np.any(np.logical_not(m.ismember(np.arange(1,md.mesh.numberofvertices+1),md.mesh.elements))):
+                       md.checkmessage("orphan nodes have been found. Check the mesh outline")
+               md = checkfield(md,'fieldname','mesh.numberofelements','>',0)
+               md = checkfield(md,'fieldname','mesh.numberofvertices','>',0)
+               md = checkfield(md,'fieldname','mesh.vertexonbase','size',[md.mesh.numberofvertices],'values',[0,1])
+               md = checkfield(md,'fieldname','mesh.vertexonsurface','size',[md.mesh.numberofvertices],'values',[0,1])
+               md = checkfield(md,'fieldname','mesh.average_vertex_connectivity','>=',9,'message',"'mesh.average_vertex_connectivity' should be at least 9 in 2d")
+
+               if solution=='ThermalSolution':
+                       md.checkmessage("thermal not supported for 2d mesh")
+
+               return md
+       # }}}
+       def domaintype(self): # {{{
+               return "2Dvertical"
+       #}}}
+       def dimension(self): # {{{
+               return 2
+       #}}}
+       def elementtype(self): # {{{
+               return "Tria"
+       #}}}
+       def vertexflags(self,value): # {{{
+               flags = np.zeros((self.numberofvertices,))
+               pos   = self.segments[np.where(self.segmentmarkers==value),0:2]-1
+               flags[pos] = 1
+               return flags
+       #}}}
+       def marshall(self,prefix,md,fid):    # {{{
+               WriteData(fid,prefix,'name','md.mesh.domain_type','data',"Domain"+self.domaintype(),'format','String');
+               WriteData(fid,prefix,'name','md.mesh.domain_dimension','data',self.dimension(),'format','Integer');
+               WriteData(fid,prefix,'name','md.mesh.elementtype','data',self.elementtype(),'format','String');
+               WriteData(fid,prefix,'object',self,'class','mesh','fieldname','x','format','DoubleMat','mattype',1)
+               WriteData(fid,prefix,'object',self,'class','mesh','fieldname','y','format','DoubleMat','mattype',1)
+               WriteData(fid,prefix,'name','md.mesh.z','data',np.zeros(self.numberofvertices),'format','DoubleMat','mattype',1);
+               WriteData(fid,prefix,'object',self,'class','mesh','fieldname','elements','format','DoubleMat','mattype',2)
+               WriteData(fid,prefix,'object',self,'class','mesh','fieldname','numberofelements','format','Integer')
+               WriteData(fid,prefix,'object',self,'class','mesh','fieldname','numberofvertices','format','Integer')
+               WriteData(fid,prefix,'object',self,'class','mesh','fieldname','vertexonbase','format','BooleanMat','mattype',1)
+               WriteData(fid,prefix,'object',self,'class','mesh','fieldname','vertexonsurface','format','BooleanMat','mattype',1)
+               WriteData(fid,prefix,'object',self,'class','mesh','fieldname','average_vertex_connectivity','format','Integer')
+       # }}}
Index: ../trunk-jpl/src/m/classes/matenhancedice.py
===================================================================
--- ../trunk-jpl/src/m/classes/matenhancedice.py        (nonexistent)
+++ ../trunk-jpl/src/m/classes/matenhancedice.py        (revision 22267)
@@ -0,0 +1,171 @@
+from fielddisplay import fielddisplay
+from project3d import project3d
+from checkfield import checkfield
+from WriteData import WriteData
+
+class matenhancedice(object):
+       """
+       MATICE class definition
+
+          Usage:
+             matenhancedice=matenhancedice();
+       """
+
+       def __init__(self): # {{{
+               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.meltingpoint              = 0.
+               self.beta                      = 0.
+               self.mixed_layer_capacity      = 0.
+               self.thermal_exchange_velocity = 0.
+               self.rheology_E                = float('NaN')
+               self.rheology_B                = float('NaN')
+               self.rheology_n                = float('NaN')
+               self.rheology_law              = ''
+
+               #giaivins: 
+               self.lithosphere_shear_modulus  = 0.
+               self.lithosphere_density        = 0.
+               self.mantle_shear_modulus       = 0.
+               self.mantle_density             = 0.
+               
+               #SLR
+               self.earth_density= 0  # average density of the Earth, (kg/m^3)
+
+               self.setdefaultparameters()
+               #}}}
+       def __repr__(self): # {{{
+               string="   Materials:"
+
+               string="%s\n%s"%(string,fielddisplay(self,"rho_ice","ice density [kg/m^3]"))
+               string="%s\n%s"%(string,fielddisplay(self,"rho_water","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_E","enhancement factor"))
+               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' or 'LliboutryDuval'"))
+               string="%s\n%s"%(string,fielddisplay(self,"lithosphere_shear_modulus","Lithosphere shear modulus [Pa]"))
+               string="%s\n%s"%(string,fielddisplay(self,"lithosphere_density","Lithosphere density [g/cm^-3]"))
+               string="%s\n%s"%(string,fielddisplay(self,"mantle_shear_modulus","Mantle shear modulus [Pa]"))
+               string="%s\n%s"%(string,fielddisplay(self,"mantle_density","Mantle density [g/cm^-3]"))
+               string="%s\n%s"%(string,fielddisplay(self,"earth_density","Mantle density [kg/m^-3]"))
+
+               return string
+               #}}}
+       def extrude(self,md): # {{{
+               self.rheology_E=project3d(md,'vector',self.rheology_E,'type','node')
+               self.rheology_B=project3d(md,'vector',self.rheology_B,'type','node')
+               self.rheology_n=project3d(md,'vector',self.rheology_n,'type','element')
+               return self
+       #}}}
+       def setdefaultparameters(self): # {{{
+               #ice density (kg/m^3)
+               self.rho_ice=917.
+
+               #ocean water density (kg/m^3)
+               self.rho_water=1023.
+
+               #fresh water density (kg/m^3)
+               self.rho_freshwater=1000.
+
+               #water viscosity (N.s/m^2)
+               self.mu_water=0.001787  
+
+               #ice heat capacity cp (J/kg/K)
+               self.heatcapacity=2093.
+
+               #ice latent heat of fusion L (J/kg)
+               self.latentheat=3.34*10**5
+
+               #ice thermal conductivity (W/m/K)
+               self.thermalconductivity=2.4
+
+               #temperate ice thermal conductivity (W/m/K)
+               self.temperateiceconductivity=0.24
+
+               #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.8*10**-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.00*10**-4
+
+               #Rheology law: what is the temperature dependence of B with T
+               #available: none, paterson and arrhenius
+               self.rheology_law='Paterson'
+
+               # GIA:
+               self.lithosphere_shear_modulus  = 6.7*10**10  # (Pa)
+               self.lithosphere_density        = 3.32        # (g/cm^-3)
+               self.mantle_shear_modulus       = 1.45*10**11 # (Pa)
+               self.mantle_density             = 3.34        # (g/cm^-3)
+               
+               #SLR
+               self.earth_density= 5512  #average density of the Earth, (kg/m^3)
+
+               return self
+               #}}}
+       def checkconsistency(self,md,solution,analyses):    # {{{
+               md = checkfield(md,'fieldname','materials.rho_ice','>',0)
+               md = checkfield(md,'fieldname','materials.rho_water','>',0)
+               md = checkfield(md,'fieldname','materials.rho_freshwater','>',0)
+               md = checkfield(md,'fieldname','materials.mu_water','>',0)
+               md = checkfield(md,'fieldname','materials.rheology_E','>',0,'timeseries',1,'NaN',1,'Inf',1)
+               md = checkfield(md,'fieldname','materials.rheology_B','>',0,'timeseries',1,'NaN',1,'Inf',1)
+               md = checkfield(md,'fieldname','materials.rheology_n','>',0,'size',[md.mesh.numberofelements])
+               md = checkfield(md,'fieldname','materials.rheology_law','values',['None','BuddJacka','Cuffey','CuffeyTemperate','Paterson','Arrhenius','LliboutryDuval'])
+
+               if 'GiaAnalysis' in analyses:
+                       md = checkfield(md,'fieldname','materials.lithosphere_shear_modulus','>',0,'numel',1)
+                       md = checkfield(md,'fieldname','materials.lithosphere_density','>',0,'numel',1)
+                       md = checkfield(md,'fieldname','materials.mantle_shear_modulus','>',0,'numel',1)
+                       md = checkfield(md,'fieldname','materials.mantle_density','>',0,'numel',1)
+               if 'SealevelriseAnalysis' in analyses:
+                       md = checkfield(md,'fieldname','materials.earth_density','>',0,'numel',1)
+               return md
+       # }}}
+       def marshall(self,prefix,md,fid):    # {{{
+               WriteData(fid,prefix,'name','md.materials.type','data',4,'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')
+               WriteData(fid,prefix,'object',self,'class','materials','fieldname','rho_freshwater','format','Double')
+               WriteData(fid,prefix,'object',self,'class','materials','fieldname','mu_water','format','Double')
+               WriteData(fid,prefix,'object',self,'class','materials','fieldname','heatcapacity','format','Double')
+               WriteData(fid,prefix,'object',self,'class','materials','fieldname','latentheat','format','Double')
+               WriteData(fid,prefix,'object',self,'class','materials','fieldname','thermalconductivity','format','Double')
+               WriteData(fid,prefix,'object',self,'class','materials','fieldname','temperateiceconductivity','format','Double')
+               WriteData(fid,prefix,'object',self,'class','materials','fieldname','meltingpoint','format','Double')
+               WriteData(fid,prefix,'object',self,'class','materials','fieldname','beta','format','Double')
+               WriteData(fid,prefix,'object',self,'class','materials','fieldname','mixed_layer_capacity','format','Double')
+               WriteData(fid,prefix,'object',self,'class','materials','fieldname','thermal_exchange_velocity','format','Double')
+               WriteData(fid,prefix,'object',self,'class','materials','fieldname','rheology_E','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts)
+               WriteData(fid,prefix,'object',self,'class','materials','fieldname','rheology_B','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts)
+               WriteData(fid,prefix,'object',self,'class','materials','fieldname','rheology_n','format','DoubleMat','mattype',2)
+               WriteData(fid,prefix,'data',self.rheology_law,'name','md.materials.rheology_law','format','String')
+
+               WriteData(fid,prefix,'object',self,'class','materials','fieldname','lithosphere_shear_modulus','format','Double')
+               WriteData(fid,prefix,'object',self,'class','materials','fieldname','lithosphere_density','format','Double','scale',10^3)
+               WriteData(fid,prefix,'object',self,'class','materials','fieldname','mantle_shear_modulus','format','Double')
+               WriteData(fid,prefix,'object',self,'class','materials','fieldname','mantle_density','format','Double','scale',10^3)
+               WriteData(fid,prefix,'object',self,'class','materials','fieldname','earth_density','format','Double')
+       # }}}
Index: ../trunk-jpl/src/m/classes/amr.py
===================================================================
--- ../trunk-jpl/src/m/classes/amr.py   (revision 22266)
+++ ../trunk-jpl/src/m/classes/amr.py   (revision 22267)
@@ -27,6 +27,7 @@
         self.deviatoricerror_resolution= 0.
         self.deviatoricerror_threshold = 0.
         self.deviatoricerror_maximum   = 0.
+       self.restart=0.
         #set defaults
         self.setdefaultparameters()
     #}}}
@@ -47,6 +48,7 @@
         string="%s\n%s"%(string,fielddisplay(self,"deviatoricerror_resolution","element length when deviatoric stress error estimator is used"))
         string="%s\n%s"%(string,fielddisplay(self,"deviatoricerror_threshold","maximum threshold deviatoricstress error permitted"))
         string="%s\n%s"%(string,fielddisplay(self,"deviatoricerror_maximum","maximum deviatoricstress error permitted"))
+       string="%s\n%s"%(string,fielddisplay(self,'restart',['indicates if ReMesh() will call before first time step']))
         return string
     #}}}
     def setdefaultparameters(self): # {{{
@@ -66,6 +68,7 @@
         self.deviatoricerror_resolution= 500.
         self.deviatoricerror_threshold = 0
         self.deviatoricerror_maximum   = 0
+       self.restart = 0.
         return self
     #}}}
     def checkconsistency(self,md,solution,analyses):    # {{{
@@ -83,6 +86,7 @@
         md = checkfield(md,'fieldname','amr.deviatoricerror_resolution','numel',[1],'>',0,'<',self.hmax,'NaN',1);
         md = checkfield(md,'fieldname','amr.deviatoricerror_threshold','numel',[1],'>=',0,'<=',1,'NaN',1);        
         md = checkfield(md,'fieldname','amr.deviatoricerror_maximum','numel',[1],'>=',0,'NaN',1,'Inf',1);
+       md = checkfield(md,'fieldname','amr.restart','numel',[1],'>=',0,'<=',1,'NaN',1)
         return md
     # }}}
     def marshall(self,prefix,md,fid):    # {{{
@@ -103,4 +107,5 @@
         WriteData(fid,prefix,'object',self,'fieldname','deviatoricerror_resolution','format','Double');
         WriteData(fid,prefix,'object',self,'fieldname','deviatoricerror_threshold','format','Double'); 
         WriteData(fid,prefix,'object',self,'fieldname','deviatoricerror_maximum','format','Double'); 
+       WriteData(fid,prefix,'object',self,'class','amr','fieldname','restart','format','Integer')
     # }}}
Index: ../trunk-jpl/src/m/classes/matestar.py
===================================================================
--- ../trunk-jpl/src/m/classes/matestar.py      (nonexistent)
+++ ../trunk-jpl/src/m/classes/matestar.py      (revision 22267)
@@ -0,0 +1,175 @@
+import numpy as np
+from fielddisplay import fielddisplay
+from project3d import project3d
+from checkfield import checkfield
+from WriteData import WriteData
+
+class matestar(object):
+       """
+       matestar class definition
+
+          Usage:
+             matestar=matestar()
+       """
+
+       def __init__(self): # {{{
+               
+               rho_ice                    = 0.
+               rho_water                  = 0.
+               rho_freshwater             = 0.
+               mu_water                   = 0.
+               heatcapacity               = 0.
+               latentheat                 = 0.
+               thermalconductivity        = 0.
+               temperateiceconductivity   = 0.
+               meltingpoint               = 0.
+               beta                       = 0.
+               mixed_layer_capacity       = 0.
+               thermal_exchange_velocity  = 0.
+               rheology_B    = float('NaN')
+               rheology_Ec   = float('NaN')
+               rheology_Es   = float('NaN')
+               rheology_law = ''
+
+               #giaivins: 
+               lithosphere_shear_modulus  = 0.
+               lithosphere_density        = 0.
+               mantle_shear_modulus       = 0.
+               mantle_density             = 0.
+
+               #slr
+               earth_density              = 0
+
+                #set default parameters:
+               self.setdefaultparameters()
+       #}}}
+       def __repr__(self): # {{{
+               string="   Materials:"
+
+               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/kg]'))
+               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/3)]'))
+               string="%s\n%s"%(string,fielddisplay(self,'rheology_Ec','compressive enhancement factor'))
+               string="%s\n%s"%(string,fielddisplay(self,'rheology_Es','shear enhancement factor'))
+               string="%s\n%s"%(string,fielddisplay(self,'rheology_law',['law for the temperature dependance of the rheology: ''None'', ''BuddJacka'', ''Cuffey'', ''CuffeyTemperate'', ''Paterson'', ''Arrhenius'' or ''LliboutryDuval''']))
+               string="%s\n%s"%(string,fielddisplay(self,'lithosphere_shear_modulus','Lithosphere shear modulus [Pa]'))
+               string="%s\n%s"%(string,fielddisplay(self,'lithosphere_density','Lithosphere density [g/cm^-3]'))
+               string="%s\n%s"%(string,fielddisplay(self,'mantle_shear_modulus','Mantle shear modulus [Pa]'))
+               string="%s\n%s"%(string,fielddisplay(self,'mantle_density','Mantle density [g/cm^-3]'))
+               string="%s\n%s"%(string,fielddisplay(self,'earth_density','Mantle density [kg/m^-3]'))
+
+               return string
+       #}}}
+       def extrude(self,md): # {{{
+               self.rheology_B=project3d(md,'vector',self.rheology_B,'type','node')
+               self.rheology_Ec=project3d(md,'vector',self.rheology_Ec,'type','node')
+               self.rheology_Es=project3d(md,'vector',self.rheology_Es,'type','node')
+                       return self
+       #}}}
+       def setdefaultparameters(self): # {{{
+               #ice density (kg/m^3)
+               self.rho_ice=917.
+
+               #ocean water density (kg/m^3)
+               self.rho_water=1023.
+
+               #fresh water density (kg/m^3)
+               self.rho_freshwater=1000.
+
+               #water viscosity (N.s/m^2)
+               self.mu_water=0.001787 
+
+               #ice heat capacity cp (J/kg/K)
+               self.heatcapacity=2093.
+
+               #ice latent heat of fusion L (J/kg)
+               self.latentheat=3.34*10**5
+
+               #ice thermal conductivity (W/m/K)
+               self.thermalconductivity=2.4
+                       
+               #wet ice thermal conductivity (W/m/K)
+               self.temperateiceconductivity=.24
+
+               #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.8*10**-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.00*10**-4
+
+               #Rheology law: what is the temperature dependence of B with T
+               #available: none, paterson and arrhenius
+               self.rheology_law='Paterson'
+
+               # GIA:
+               self.lithosphere_shear_modulus  = 6.7*10**10  # (Pa)
+               self.lithosphere_density        = 3.32      # (g/cm^-3)
+               self.mantle_shear_modulus       = 1.45*10**11 # (Pa)
+               self.mantle_density             = 3.34      # (g/cm^-3)
+
+               #SLR
+               self.earth_density= 5512  # average density of the Earth, (kg/m^3)
+
+               return self
+       #}}}
+       def checkconsistency(self,md,solution,analyses):    # {{{
+               md = checkfield(md,'fieldname','materials.rho_ice','>',0)
+               md = checkfield(md,'fieldname','materials.rho_water','>',0)
+               md = checkfield(md,'fieldname','materials.rho_freshwater','>',0)
+               md = checkfield(md,'fieldname','materials.mu_water','>',0)
+               md = checkfield(md,'fieldname','materials.rheology_B','>',0,'size',[md.mesh.numberofvertices],'NaN',1,'Inf',1)
+               md = checkfield(md,'fieldname','materials.rheology_Ec','>',0,'size',[md.mesh.numberofvertices],'NaN',1,'Inf',1)
+               md = checkfield(md,'fieldname','materials.rheology_Es','>',0,'size',[md.mesh.numberofvertices],'NaN',1,'Inf',1)
+               md = checkfield(md,'fieldname','materials.rheology_law','values',['None','BuddJacka', 'Cuffey','CuffeyTemperate','Paterson','Arrhenius','LliboutryDuval'])
+
+               if 'GiaAnalysis' in analyses:
+                       md = checkfield(md,'fieldname','materials.lithosphere_shear_modulus','>',0,'numel',1)
+                       md = checkfield(md,'fieldname','materials.lithosphere_density','>',0,'numel',1)
+                       md = checkfield(md,'fieldname','materials.mantle_shear_modulus','>',0,'numel',1)
+                       md = checkfield(md,'fieldname','materials.mantle_density','>',0,'numel',1)
+               if 'SealevelriseAnalysis' in analyses:
+                       md = checkfield(md,'fieldname','materials.earth_density','>',0,'numel',1)
+
+               return md
+       # }}}
+       def marshall(self,prefix,md,fid):    # {{{
+               WriteData(fid,prefix,'name','md.materials.type','data',2,'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')
+               WriteData(fid,prefix,'object',self,'class','materials','fieldname','rho_freshwater','format','Double')
+               WriteData(fid,prefix,'object',self,'class','materials','fieldname','mu_water','format','Double')
+               WriteData(fid,prefix,'object',self,'class','materials','fieldname','heatcapacity','format','Double')
+               WriteData(fid,prefix,'object',self,'class','materials','fieldname','latentheat','format','Double')
+               WriteData(fid,prefix,'object',self,'class','materials','fieldname','thermalconductivity','format','Double')
+               WriteData(fid,prefix,'object',self,'class','materials','fieldname','temperateiceconductivity','format','Double')
+               WriteData(fid,prefix,'object',self,'class','materials','fieldname','meltingpoint','format','Double')
+               WriteData(fid,prefix,'object',self,'class','materials','fieldname','beta','format','Double')
+               WriteData(fid,prefix,'object',self,'class','materials','fieldname','mixed_layer_capacity','format','Double')
+               WriteData(fid,prefix,'object',self,'class','materials','fieldname','thermal_exchange_velocity','format','Double')
+               WriteData(fid,prefix,'object',self,'class','materials','fieldname','rheology_B','format','DoubleMat','mattype',1)
+               WriteData(fid,prefix,'object',self,'class','materials','fieldname','rheology_Ec','format','DoubleMat','mattype',1)
+               WriteData(fid,prefix,'object',self,'class','materials','fieldname','rheology_Es','format','DoubleMat','mattype',1)
+               WriteData(fid,prefix,'data',self.rheology_law,'name','md.materials.rheology_law','format','String')
+
+               WriteData(fid,prefix,'object',self,'class','materials','fieldname','lithosphere_shear_modulus','format','Double')
+               WriteData(fid,prefix,'object',self,'class','materials','fieldname','lithosphere_density','format','Double','scale',10^3)
+               WriteData(fid,prefix,'object',self,'class','materials','fieldname','mantle_shear_modulus','format','Double')
+               WriteData(fid,prefix,'object',self,'class','materials','fieldname','mantle_density','format','Double','scale',10**3)
+               WriteData(fid,prefix,'object',self,'class','materials','fieldname','earth_density','format','Double')
+       # }}}
Index: ../trunk-jpl/src/m/classes/frictionsommers.py
===================================================================
--- ../trunk-jpl/src/m/classes/frictionsommers.py       (nonexistent)
+++ ../trunk-jpl/src/m/classes/frictionsommers.py       (revision 22267)
@@ -0,0 +1,47 @@
+from fielddisplay import fielddisplay
+from project3d import project3d
+from checkfield import checkfield
+from WriteData import WriteData
+
+class frictionsommers(object):
+    """
+    FRICTIONSOMMERS class definition
+
+    Usage:
+        friction=frictionsommers()
+    """
+
+    def __init__(self,md): # {{{
+        self.coefficient = md.friction.coefficient
+       #set defaults
+       self.setdefaultparameters()
+
+    #}}}
+    def __repr__(self): # {{{
+       string="Basal shear stress parameters: Sigma_b = coefficient^2 * Neff * u_b\n(effective stress Neff=rho_ice*g*thickness+rho_water*g*(head-b))"
+
+       string="%s\n%s"%(string,fielddisplay(self,"coefficient","friction coefficient [SI]"))
+       return string
+    #}}}
+    def extrude(self,md): # {{{
+       self.coefficient=project3d(md,'vector',self.coefficient,'type','node','layer',1)        
+       return self
+    #}}}
+    def setdefaultparameters(self): # {{{
+       return self
+    #}}}
+    def checkconsistency(self,md,solution,analyses):    # {{{
+
+       #Early return
+       if 'StressbalanceAnalysis' not in analyses and 'ThermalAnalysis' not in analyses:
+           return md
+
+       md = checkfield(md,'fieldname','friction.coefficient','timeseries',1,'NaN',1,'Inf',1)
+       return md
+
+    # }}}
+    def marshall(self,prefix,md,fid):    # {{{
+       yts=md.constants.yts
+       WriteData(fid,prefix,'name','md.friction.law','data',8,'format','Integer')
+       WriteData(fid,prefix,'object',self,'fieldname','coefficient','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts)     
+    # }}}
Index: ../trunk-jpl/src/m/classes/SMBgemb.py
===================================================================
--- ../trunk-jpl/src/m/classes/SMBgemb.py       (nonexistent)
+++ ../trunk-jpl/src/m/classes/SMBgemb.py       (revision 22267)
@@ -0,0 +1,368 @@
+import numpy as np
+from fielddisplay import fielddisplay
+from checkfield import checkfield
+from WriteData import WriteData
+from project3d import project3d
+
+class SMBgemb(object):
+       """
+       SMBgemb Class definition
+
+          Usage:
+             SMB = SMBgemb()
+       """
+
+       def __init__(self): # {{{
+               #each one of these properties is a transient forcing to the GEMB model, loaded from meteorological data derived 
+               #from an automatic weather stations (AWS). Each property is therefore a matrix, of size (numberofvertices x number 
+               #of time steps. )
+
+               #solution choices
+               #check these:
+               #isgraingrowth
+               #isalbedo
+               #isshortwave
+               #isthermal
+               #isaccumulation
+               #ismelt
+               #isdensification
+               #isturbulentflux    
+
+               #inputs: 
+               Ta    = float('NaN')    #2 m air temperature, in Kelvin
+               V     = float('NaN')    #wind speed (m/s-1)
+               dswrf = float('NaN')    #downward shortwave radiation flux [W/m^2]
+               dlwrf = float('NaN')    #downward longwave radiation flux [W/m^2]
+               P     = float('NaN')    #precipitation [mm w.e. / m^2]
+               eAir  = float('NaN')    #screen level vapor pressure [Pa]
+               pAir  = float('NaN')    #surface pressure [Pa]
+               
+               Tmean = float('NaN')    #mean annual temperature [K]
+               C     = float('NaN')    #mean annual snow accumulation [kg m-2 yr-1]
+               Tz    = float('NaN')    #height above ground at which temperature (T) was sampled [m]
+               Vz    = float('NaN')    #height above ground at which wind (V) eas sampled [m]
+        
+               # Initialization of snow properties
+               Dzini = float('NaN')    #cell depth (m)
+               Dini = float('NaN')     #snow density (kg m-3)
+               Reini = float('NaN')    #effective grain size (mm)
+               Gdnini = float('NaN')   #grain dricity (0-1)
+               Gspini = float('NaN')   #grain sphericity (0-1)
+               ECini = float('NaN')    #evaporation/condensation (kg m-2)
+               Wini = float('NaN')     #Water content (kg m-2)
+               Aini = float('NaN')     #albedo (0-1)
+               Tini = float('NaN')     #snow temperature (K)
+               Sizeini = float('NaN')  #Number of layers
+
+               #settings: 
+               aIdx   = float('NaN')   #method for calculating albedo and subsurface absorption (default is 1)
+                                       # 1: effective grain radius [Gardner & Sharp, 2009]
+                                         # 2: effective grain radius [Brun et al., 2009]
+                                         # 3: density and cloud amount [Greuell & Konzelmann, 1994]
+                                         # 4: exponential time decay & wetness [Bougamont & Bamber, 2005]
+               swIdx  = float('NaN')   # apply all SW to top grid cell (0) or allow SW to penetrate surface (1) (default 1)
+
+               denIdx = float('NaN')   #densification model to use (default is 2):
+                                       # 1 = emperical model of Herron and Langway (1980)
+                                       # 2 = semi-emerical model of Anthern et al. (2010)
+                                       # 3 = DO NOT USE: physical model from Appix B of Anthern et al. (2010)
+                                       # 4 = DO NOT USE: emperical model of Li and Zwally (2004)
+                                       # 5 = DO NOT USE: modified emperical model (4) by Helsen et al. (2008)
+
+               zTop  = float('NaN')    # depth over which grid length is constant at the top of the snopack (default 10) [m]
+               dzTop = float('NaN')    # initial top vertical grid spacing (default .05) [m] 
+               dzMin = float('NaN')    # initial min vertical allowable grid spacing (default dzMin/2) [m] 
+               zY    = float('NaN')    # strech grid cells bellow top_z by a [top_dz * y ^ (cells bellow top_z)]
+               zMax = float('NaN')     #initial max model depth (default is min(thickness,500)) [m]
+               zMin = float('NaN')     #initial min model depth (default is min(thickness,30)) [m]
+               outputFreq = float('NaN')       #output frequency in days (default is monthly, 30)
+
+               #specific albedo parameters: 
+               #Method 1 and 2: 
+               aSnow = float('NaN')    # new snow albedo (0.64 - 0.89)
+               aIce  = float('NaN')    # range 0.27-0.58 for old snow
+                       #Method 3: Radiation Correction Factors -> only used for met station data and Greuell & Konzelmann, 1994 albedo
+               cldFrac = float('NaN')  # average cloud amount
+                       #Method 4: additonal tuning parameters albedo as a funtion of age and water content (Bougamont et al., 2005)
+               t0wet = float('NaN')    # time scale for wet snow (15-21.9) 
+               t0dry = float('NaN')    # warm snow timescale (30) 
+               K     = float('NaN')    # time scale temperature coef. (7) 
+
+               #densities:
+               InitDensityScaling =  float('NaN')      #initial scaling factor multiplying the density of ice, which describes the density of the snowpack.
+               
+               requested_outputs      = []
+
+               #Several fields are missing from the standard GEMB model, which are capture intrinsically by ISSM. 
+               #dateN: that's the last row of the above fields. 
+               #dt:    included in dateN. Not an input.  
+               #elev:  this is taken from the ISSM surface itself.
+
+               #}}}
+       def __repr__(self): # {{{
+               #string = "   surface forcings parameters:"
+               #string = "#s\n#s"%(string,fielddisplay(self,'mass_balance','surface mass balance [m/yr ice eq]'))
+               #string = "#s\n#s"%(string,fielddisplay(self,'requested_outputs','additional outputs requested'))
+               string = '   surface forcings for SMB GEMB model :'
+                       
+               string="%s\n%s"%(string,fielddisplay(self,'issmbgradients','is smb gradients method activated (0 or 1, default is 0)'))
+               string = "%s\n%s"%(string,fielddisplay(self,'isgraingrowth','run grain growth module (default true)'))
+               string = "%s\n%s"%(string,fielddisplay(self,'isalbedo','run albedo module (default true)'))
+               string = "%s\n%s"%(string,fielddisplay(self,'isshortwave','run short wave module (default true)'))
+               string = "%s\n%s"%(string,fielddisplay(self,'isthermal','run thermal module (default true)'))
+               string = "%s\n%s"%(string,fielddisplay(self,'isaccumulation','run accumulation module (default true)'))
+               string = "%s\n%s"%(string,fielddisplay(self,'ismelt','run melting  module (default true)'))
+               string = "%s\n%s"%(string,fielddisplay(self,'isdensification','run densification module (default true)'))
+               string = "%s\n%s"%(string,fielddisplay(self,'isturbulentflux','run turbulant heat fluxes module (default true)'))
+               string = "%s\n%s"%(string,fielddisplay(self,'Ta','2 m air temperature, in Kelvin'))
+               string = "%s\n%s"%(string,fielddisplay(self,'V','wind speed (m/s-1)'))
+               string = "%s\n%s"%(string,fielddisplay(self,'dlwrf','downward shortwave radiation flux [W/m^2]'))
+               string = "%s\n%s"%(string,fielddisplay(self,'dswrf','downward longwave radiation flux [W/m^2]'))
+               string = "%s\n%s"%(string,fielddisplay(self,'P','precipitation [mm w.e. / m^2]'))
+               string = "%s\n%s"%(string,fielddisplay(self,'eAir','screen level vapor pressure [Pa]'))
+               string = "%s\n%s"%(string,fielddisplay(self,'pAir','surface pressure [Pa]'))
+               string = "%s\n%s"%(string,fielddisplay(self,'Tmean','mean annual temperature [K]'))
+               string = "%s\n%s"%(string,fielddisplay(self,'C','mean annual snow accumulation [kg m-2 yr-1]'))
+               string = "%s\n%s"%(string,fielddisplay(self,'Tz','height above ground at which temperature (T) was sampled [m]'))
+               string = "%s\n%s"%(string,fielddisplay(self,'Vz','height above ground at which wind (V) eas sampled [m]'))
+               string = "%s\n%s"%(string,fielddisplay(self,'zTop','depth over which grid length is constant at the top of the snopack (default 10) [m]'))
+               string = "%s\n%s"%(string,fielddisplay(self,'dzTop','initial top vertical grid spacing (default .05) [m] '))
+               string = "%s\n%s"%(string,fielddisplay(self,'dzMin','initial min vertical allowable grid spacing (default dzMin/2) [m] '))
+               string = "%s\n%s"%(string,fielddisplay(self,'zMax','initial max model depth (default is min(thickness,500)) [m]'))
+               string = "%s\n%s"%(string,fielddisplay(self,'zMin','initial min model depth (default is min(thickness,30)) [m]'))
+               string = "%s\n%s"%(string,fielddisplay(self,'zY','strech grid cells bellow top_z by a [top_dz * y ^ (cells bellow top_z)]'))
+               string = "%s\n%s"%(string,fielddisplay(self,'InitDensityScaling',['initial scaling factor multiplying the density of ice','which describes the density of the snowpack.']))
+               string = "%s\n%s"%(string,fielddisplay(self,'outputFreq','output frequency in days (default is monthly, 30)'))
+               string = "%s\n%s"%(string,fielddisplay(self,'aIdx',['method for calculating albedo and subsurface absorption (default is 1)',
+                                               '1: effective grain radius [Gardner & Sharp, 2009]',
+                                               '2: effective grain radius [Brun et al., 2009]',
+                                               '3: density and cloud amount [Greuell & Konzelmann, 1994]',
+                                               '4: exponential time decay & wetness [Bougamont & Bamber, 2005]']))
+                               
+               #snow properties init
+               string = "%s\n%s"%(string,fielddisplay(self,'Dzini','Initial cel depth when restart [m]'))
+               string = "%s\n%s"%(string,fielddisplay(self,'Dini','Initial snow density when restart [kg m-3]'))
+               string = "%s\n%s"%(string,fielddisplay(self,'Reini','Initial grain size when restart [mm]'))
+               string = "%s\n%s"%(string,fielddisplay(self,'Gdnini','Initial grain dricity when restart [-]'))
+               string = "%s\n%s"%(string,fielddisplay(self,'Gspini','Initial grain sphericity when restart [-]'))
+               string = "%s\n%s"%(string,fielddisplay(self,'ECini','Initial evaporation/condensation when restart [kg m-2]'))
+               string = "%s\n%s"%(string,fielddisplay(self,'Wini','Initial snow water content when restart [kg m-2]'))
+               string = "%s\n%s"%(string,fielddisplay(self,'Aini','Initial albedo when restart [-]'))
+               string = "%s\n%s"%(string,fielddisplay(self,'Tini','Initial snow temperature when restart [K]'))
+               string = "%s\n%s"%(string,fielddisplay(self,'Sizeini','Initial number of layers when restart [K]'))
+                        
+               #additional albedo parameters: 
+               if type(self.aIdx) == list or type(self.aIdx) == type(np.array([1,2])) and (self.aIdx == [1,2] or (1 in self.aIdx and 2 in self.aIdx)):
+                       string = "%s\n%s"%(string,fielddisplay(self,'aSnow','new snow albedo (0.64 - 0.89)'))
+                       string = "%s\n%s"%(string,fielddisplay(self,'aIce','albedo of ice (0.27-0.58)'))
+               elif self.aIdx == 3:
+                       string = "%s\n%s"%(string,fielddisplay(self,'cldFrac','average cloud amount'))
+               elif self.aIdx == 4:
+                       string = "%s\n%s"%(string,fielddisplay(self,'t0wet','time scale for wet snow (15-21.9) [d]'))
+                       string = "%s\n%s"%(string,fielddisplay(self,'t0dry','warm snow timescale (30) [d]'))
+                       string = "%s\n%s"%(string,fielddisplay(self,'K','time scale temperature coef. (7) [d]'))
+               
+               string = "%s\n%s"%(string,fielddisplay(self,'swIdx','apply all SW to top grid cell (0) or allow SW to penetrate surface (1) [default 1]'))
+               string = "%s\n%s"%(string,fielddisplay(self,'denIdx',['densification model to use (default is 2):',
+                                               '1 = emperical model of Herron and Langway (1980)',
+                                               '2 = semi-emerical model of Anthern et al. (2010)',
+                                               '3 = DO NOT USE: physical model from Appix B of Anthern et al. (2010)',
+                                               '4 = DO NOT USE: emperical model of Li and Zwally (2004)',
+                                               '5 = DO NOT USE: modified emperical model (4) by Helsen et al. (2008)']))
+               string = "%s\n%s"%(string,fielddisplay(self,'requested_outputs','additional outputs requested'))
+               return string
+       #}}}
+
+       def extrude(self,md): # {{{
+
+               self.Ta = project3d(md,'vector',self.Ta,'type','node')
+               self.V = project3d(md,'vector',self.V,'type','node')
+               self.dswrf = project3d(md,'vector',self.dswrf,'type','node')
+               self.dswrf = project3d(md,'vector',self.dswrf,'type','node')
+               self.P = project3d(md,'vector',self.P,'type','node')
+               self.eAir = project3d(md,'vector',self.eAir,'type','node')
+               self.pAir = project3d(md,'vector',self.pAir,'type','node')
+               return self
+       #}}}
+       def defaultoutputs(self,md): # {{{
+               return ['SmbMassBalance']
+       #}}}
+
+       def setdefaultparameters(self,mesh,geometry): # {{{
+               self.isgraingrowth = 1
+               self.isalbedo = 1
+               self.isshortwave = 1
+               self.isthermal = 1
+               self.isaccumulation = 1
+               self.ismelt = 1
+               self.isdensification = 1
+               self.isturbulentflux = 1
+       
+               self.aIdx = 1
+               self.swIdx = 1
+               self.denIdx = 2
+               self.zTop = 10*np.ones((mesh.numberofelements,))
+               self.dzTop = .05* np.ones((mesh.numberofelements,))
+               self.dzMin = self.dzTop/2
+               self.InitDensityScaling = 1.0
+               
+               self.zMax = 250*np.ones((mesh.numberofelements,))
+               self.zMin = 130*np.ones((mesh.numberofelements,))
+               self.zY = 1.10*np.ones((mesh.numberofelements,))
+               self.outputFreq = 30
+               
+               #additional albedo parameters
+               self.aSnow = 0.85
+               self.aIce = 0.48
+               self.cldFrac = 0.1 
+               self.t0wet = 15
+               self.t0dry = 30
+               self.K = 7
+        
+               self.Dzini = 0.05*np.ones((mesh.numberofelements,2))
+               self.Dini = 910.0*np.ones((mesh.numberofelements,2)) 
+               self.Reini = 2.5*np.ones((mesh.numberofelements,2))
+               self.Gdnini = 0.0*np.ones((mesh.numberofelements,2))
+               self.Gspini = 0.0*np.ones((mesh.numberofelements,2))
+               self.ECini = 0.0*np.ones((mesh.numberofelements,))
+               self.Wini = 0.0*np.ones((mesh.numberofelements,2))
+               self.Aini = self.aSnow*np.ones((mesh.numberofelements,2))
+               self.Tini = 273.15*np.ones((mesh.numberofelements,2))
+#              /!\ Default value of Tini must be equal to Tmean but don't know Tmean yet (computed when atmospheric forcings are interpolated on mesh). 
+#              If initialization without restart, this value will be overwritten when snow parameters are retrieved in Element.cpp 
+               self.Sizeini = 2*np.ones((mesh.numberofelements,))
+       #}}}
+
+       def checkconsistency(self,md,solution,analyses):    # {{{
+
+               md = checkfield(md,'fieldname','smb.isgraingrowth','values',[0,1])
+               md = checkfield(md,'fieldname','smb.isalbedo','values',[0,1])
+               md = checkfield(md,'fieldname','smb.isshortwave','values',[0,1])
+               md = checkfield(md,'fieldname','smb.isthermal','values',[0,1])
+               md = checkfield(md,'fieldname','smb.isaccumulation','values',[0,1])
+               md = checkfield(md,'fieldname','smb.ismelt','values',[0,1])
+               md = checkfield(md,'fieldname','smb.isdensification','values',[0,1])
+               md = checkfield(md,'fieldname','smb.isturbulentflux','values',[0,1])
+
+               md = checkfield(md,'fieldname','smb.Ta','timeseries',1,'NaN',1,'Inf',1,'>',273-100,'<',273+100) #-100/100 celsius min/max value
+               md = checkfield(md,'fieldname','smb.V','timeseries',1,'NaN',1,'Inf',1,'> = ',0,'<',45) #max 500 km/h
+               md = checkfield(md,'fieldname','smb.dswrf','timeseries',1,'NaN',1,'Inf',1,'> = ',0,'< = ',1400)
+               md = checkfield(md,'fieldname','smb.dlwrf','timeseries',1,'NaN',1,'Inf',1,'> = ',0)
+               md = checkfield(md,'fieldname','smb.P','timeseries',1,'NaN',1,'Inf',1,'> = ',0,'< = ',100)
+               md = checkfield(md,'fieldname','smb.eAir','timeseries',1,'NaN',1,'Inf',1)
+
+               md = checkfield(md,'fieldname','smb.Tmean','size',[md.mesh.numberofelements],'NaN',1,'Inf',1,'>',273-100,'<',273+100) #-100/100 celsius min/max value
+               md = checkfield(md,'fieldname','smb.C','size',[md.mesh.numberofelements],'NaN',1,'Inf',1,'> = ',0) 
+               md = checkfield(md,'fieldname','smb.Tz','size',[md.mesh.numberofelements],'NaN',1,'Inf',1,'> = ',0,'< = ',5000) 
+               md = checkfield(md,'fieldname','smb.Vz','size',[md.mesh.numberofelements],'NaN',1,'Inf',1,'> = ',0,'< = ',5000) 
+
+               md = checkfield(md,'fieldname','smb.aIdx','NaN',1,'Inf',1,'values',[1,2,3,4])
+               md = checkfield(md,'fieldname','smb.swIdx','NaN',1,'Inf',1,'values',[0,1])
+               md = checkfield(md,'fieldname','smb.denIdx','NaN',1,'Inf',1,'values',[1,2,3,4,5])
+
+               md = checkfield(md,'fieldname','smb.zTop','NaN',1,'Inf',1,'> = ',0)
+               md = checkfield(md,'fieldname','smb.dzTop','NaN',1,'Inf',1,'>',0)
+               md = checkfield(md,'fieldname','smb.dzMin','NaN',1,'Inf',1,'>',0)
+               md = checkfield(md,'fieldname','smb.zY','NaN',1,'Inf',1,'> = ',1)
+               md = checkfield(md,'fieldname','smb.outputFreq','NaN',1,'Inf',1,'>',0,'<',10*365) #10 years max 
+               md = checkfield(md,'fieldname','smb.InitDensityScaling','NaN',1,'Inf',1,'> = ',0,'< = ',1)
+
+               if type(self.aIdx) == list or type(self.aIdx) == type(np.array([1,2])) and (self.aIdx == [1,2] or (1 in self.aIdx and 2 in self.aIdx)):
+                       md = checkfield(md,'fieldname','smb.aSnow','NaN',1,'Inf',1,'> = ',.64,'< = ',.89)
+                       md = checkfield(md,'fieldname','smb.aIce','NaN',1,'Inf',1,'> = ',.27,'< = ',.58)
+               elif self.aIdx == 3:
+                       md = checkfield(md,'fieldname','smb.cldFrac','NaN',1,'Inf',1,'> = ',0,'< = ',1)
+               elif self.aIdx == 4:
+                       md = checkfield(md,'fieldname','smb.t0wet','NaN',1,'Inf',1,'> = ',15,'< = ',21.9)
+                       md = checkfield(md,'fieldname','smb.t0dry','NaN',1,'Inf',1,'> = ',30,'< = ',30)
+                       md = checkfield(md,'fieldname','smb.K','NaN',1,'Inf',1,'> = ',7,'< = ',7)
+                       
+
+               #check zTop is < local thickness:
+               he = np.sum(md.geometry.thickness[md.mesh.elements-1],axis=1)/np.size(md.mesh.elements,1)
+               if np.any(he<self.zTop):
+                       error('SMBgemb consistency check error: zTop should be smaller than local ice thickness')
+               
+               md = checkfield(md,'fieldname','smb.requested_outputs','stringrow',1)
+               return md
+       # }}}
+
+       def marshall(self,prefix,md,fid):    # {{{
+
+               yts = md.constants.yts
+
+               WriteData(fid,prefix,'name','md.smb.model','data',8,'format','Integer')
+                       
+               WriteData(fid,prefix,'object',self,'class','smb','fieldname','isgraingrowth','format','Boolean')
+               WriteData(fid,prefix,'object',self,'class','smb','fieldname','isalbedo','format','Boolean')
+               WriteData(fid,prefix,'object',self,'class','smb','fieldname','isshortwave','format','Boolean')
+               WriteData(fid,prefix,'object',self,'class','smb','fieldname','isthermal','format','Boolean')
+               WriteData(fid,prefix,'object',self,'class','smb','fieldname','isaccumulation','format','Boolean')
+               WriteData(fid,prefix,'object',self,'class','smb','fieldname','ismelt','format','Boolean')
+               WriteData(fid,prefix,'object',self,'class','smb','fieldname','isdensification','format','Boolean')
+               WriteData(fid,prefix,'object',self,'class','smb','fieldname','isturbulentflux','format','Boolean')
+            
+               WriteData(fid,prefix,'object',self,'class','smb','fieldname','Ta','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts)
+               WriteData(fid,prefix,'object',self,'class','smb','fieldname','V','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts)
+               WriteData(fid,prefix,'object',self,'class','smb','fieldname','dswrf','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts)
+               WriteData(fid,prefix,'object',self,'class','smb','fieldname','dlwrf','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts)
+               WriteData(fid,prefix,'object',self,'class','smb','fieldname','P','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts)
+               WriteData(fid,prefix,'object',self,'class','smb','fieldname','eAir','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts)
+               WriteData(fid,prefix,'object',self,'class','smb','fieldname','pAir','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts)         
+
+               WriteData(fid,prefix,'object',self,'class','smb','fieldname','Tmean','format','DoubleMat','mattype',2)
+               WriteData(fid,prefix,'object',self,'class','smb','fieldname','C','format','DoubleMat','mattype',2)
+               WriteData(fid,prefix,'object',self,'class','smb','fieldname','Tz','format','DoubleMat','mattype',2)
+               WriteData(fid,prefix,'object',self,'class','smb','fieldname','Vz','format','DoubleMat','mattype',2)
+               WriteData(fid,prefix,'object',self,'class','smb','fieldname','zTop','format','DoubleMat','mattype',2)
+               WriteData(fid,prefix,'object',self,'class','smb','fieldname','dzTop','format','DoubleMat','mattype',2)
+               WriteData(fid,prefix,'object',self,'class','smb','fieldname','dzMin','format','DoubleMat','mattype',2)
+               WriteData(fid,prefix,'object',self,'class','smb','fieldname','zY','format','DoubleMat','mattype',2)
+               WriteData(fid,prefix,'object',self,'class','smb','fieldname','zMax','format','DoubleMat','mattype',2)
+               WriteData(fid,prefix,'object',self,'class','smb','fieldname','zMin','format','DoubleMat','mattype',2)
+               
+               WriteData(fid,prefix,'object',self,'class','smb','fieldname','aIdx','format','Integer')
+               WriteData(fid,prefix,'object',self,'class','smb','fieldname','swIdx','format','Integer')
+               WriteData(fid,prefix,'object',self,'class','smb','fieldname','denIdx','format','Integer')
+               WriteData(fid,prefix,'object',self,'class','smb','fieldname','InitDensityScaling','format','Double')
+               WriteData(fid,prefix,'object',self,'class','smb','fieldname','outputFreq','format','Double')
+               WriteData(fid,prefix,'object',self,'class','smb','fieldname','aSnow','format','Double')
+               WriteData(fid,prefix,'object',self,'class','smb','fieldname','aIce','format','Double')
+               WriteData(fid,prefix,'object',self,'class','smb','fieldname','cldFrac','format','Double')
+               WriteData(fid,prefix,'object',self,'class','smb','fieldname','t0wet','format','Double')
+               WriteData(fid,prefix,'object',self,'class','smb','fieldname','t0dry','format','Double')
+               WriteData(fid,prefix,'object',self,'class','smb','fieldname','K','format','Double')
+            
+               #snow properties init
+               WriteData(fid,prefix,'object',self,'class','smb','fieldname','Dzini','format','DoubleMat','mattype',3)
+               WriteData(fid,prefix,'object',self,'class','smb','fieldname','Dini','format','DoubleMat','mattype',3)
+               WriteData(fid,prefix,'object',self,'class','smb','fieldname','Reini','format','DoubleMat','mattype',3)
+               WriteData(fid,prefix,'object',self,'class','smb','fieldname','Gdnini','format','DoubleMat','mattype',3)
+               WriteData(fid,prefix,'object',self,'class','smb','fieldname','Gspini','format','DoubleMat','mattype',3)
+               WriteData(fid,prefix,'object',self,'class','smb','fieldname','ECini','format','DoubleMat','mattype',2)
+               WriteData(fid,prefix,'object',self,'class','smb','fieldname','Wini','format','DoubleMat','mattype',3)
+               WriteData(fid,prefix,'object',self,'class','smb','fieldname','Aini','format','DoubleMat','mattype',3)
+               WriteData(fid,prefix,'object',self,'class','smb','fieldname','Tini','format','DoubleMat','mattype',3)
+               WriteData(fid,prefix,'object',self,'class','smb','fieldname','Sizeini','format','IntMat','mattype',2)
+
+               #figure out dt from forcings: 
+               time = self.Ta[-1] #assume all forcings are on the same time step
+               dtime = np.diff(time,n=1,axis=0)
+               dt = min(dtime) / yts
+            
+               WriteData(fid,prefix,'data',dt,'name','md.smb.dt','format','Double','scale',yts)
+
+               # Check if smb_dt goes evenly into transient core time step
+               if (md.timestepping.time_step % dt >=  1e-10):
+                       error('smb_dt/dt = #f. The number of SMB time steps in one transient core time step has to be an an integer',md.timestepping.time_step/dt)
+                       
+               #process requested outputs
+               outputs = self.requested_outputs
+               indices = [i for i, x in enumerate(outputs) if x == 'default']
+               if len(indices) > 0:
+                       outputscopy=outputs[0:max(0,indices[0]-1)]+self.defaultoutputs(md)+outputs[indices[0]+1:]
+                       outputs    =outputscopy
+               
+               WriteData(fid,prefix,'data',outputs,'name','md.smb.requested_outputs','format','StringArray')
+       # }}}
+
Index: ../trunk-jpl/src/m/classes/taoinversion.py
===================================================================
--- ../trunk-jpl/src/m/classes/taoinversion.py  (revision 22266)
+++ ../trunk-jpl/src/m/classes/taoinversion.py  (revision 22267)
@@ -5,31 +5,33 @@
 from fielddisplay import fielddisplay
 from IssmConfig import IssmConfig
 from marshallcostfunctions import marshallcostfunctions
+from supportedcontrols import *
+from supportedcostfunctions import *
 
-
-class taoinversion:
+class taoinversion(object):
        def __init__(self):
-               iscontrol                   = 0
-               incomplete_adjoint          = 0
-               control_parameters          = float('NaN')
-               maxsteps                    = 0
-               maxiter                     = 0
-               fatol                       = 0
-               frtol                       = 0
-               gatol                       = 0
-               grtol                       = 0
-               gttol                       = 0
-               algorithm                   = ''
-               cost_functions              = float('NaN')
-               cost_functions_coefficients = float('NaN')
-               min_parameters              = float('NaN')
-               max_parameters              = float('NaN')
-               vx_obs                      = float('NaN')
-               vy_obs                      = float('NaN')
-               vz_obs                      = float('NaN')
-               vel_obs                     = float('NaN')
-               thickness_obs               = float('NaN')
-               surface_obs                 = float('NaN')
+               self.iscontrol                   = 0
+               self.incomplete_adjoint          = 0
+               self.control_parameters          = float('NaN')
+               self.maxsteps                    = 0
+               self.maxiter                     = 0
+               self.fatol                       = 0
+               self.frtol                       = 0
+               self.gatol                       = 0
+               self.grtol                       = 0
+               self.gttol                       = 0
+               self.algorithm                   = ''
+               self.cost_functions              = float('NaN')
+               self.cost_functions_coefficients = float('NaN')
+               self.min_parameters              = float('NaN')
+               self.max_parameters              = float('NaN')
+               self.vx_obs                      = float('NaN')
+               self.vy_obs                      = float('NaN')
+               self.vz_obs                      = float('NaN')
+               self.vel_obs                     = float('NaN')
+               self.thickness_obs               = float('NaN')
+               self.surface_obs                 = float('NaN')
+               self.setdefaultparameters()
 
        def __repr__(self):
                string = '   taoinversion parameters:'
@@ -97,7 +99,6 @@
                
                #several responses can be used:
                self.cost_functions=101;
-
                return self
 
        def extrude(self,md):
@@ -118,17 +119,17 @@
                return self
 
        def checkconsistency(self,md,solution,analyses):
-               if not self.control:
+               if not self.iscontrol:
                        return md
                if not IssmConfig('_HAVE_TAO_')[0]:
                        md = checkmessage(md,['TAO has not been installed, ISSM needs to be reconfigured and recompiled with TAO'])
 
 
-               num_controls= np.numel(md.inversion.control_parameters)
-               num_costfunc= np.size(md.inversion.cost_functions,2)
+               num_controls= np.size(md.inversion.control_parameters)
+               num_costfunc= np.size(md.inversion.cost_functions)
 
                md = checkfield(md,'fieldname','inversion.iscontrol','values',[0, 1])
-               md = checkfield(md,'fieldname','inversion.incomplete_adjoint','values',[0, 1])
+               md = checkfield(md,'fieldname','inversion.incomplete_adjoint','values',[0,1])
                md = checkfield(md,'fieldname','inversion.control_parameters','cell',1,'values',supportedcontrols())
                md = checkfield(md,'fieldname','inversion.maxsteps','numel',1,'>=',0)
                md = checkfield(md,'fieldname','inversion.maxiter','numel',1,'>=',0)
@@ -142,12 +143,12 @@
                PETSCMAJOR = IssmConfig('_PETSC_MAJOR_')[0]
                PETSCMINOR = IssmConfig('_PETSC_MINOR_')[0]
                if(PETSCMAJOR>3 or (PETSCMAJOR==3 and PETSCMINOR>=5)):
-                       md = checkfield(md,'fieldname','inversion.algorithm','values',{'blmvm','cg','lmvm'})
+                       md = checkfield(md,'fieldname','inversion.algorithm','values',['blmvm','cg','lmvm'])
                else:
-                       md = checkfield(md,'fieldname','inversion.algorithm','values',{'tao_blmvm','tao_cg','tao_lmvm'})
+                       md = checkfield(md,'fieldname','inversion.algorithm','values',['tao_blmvm','tao_cg','tao_lmvm'])
 
 
-               md = checkfield(md,'fieldname','inversion.cost_functions','size',[1, num_costfunc],'values',supportedcostfunctions())
+               md = checkfield(md,'fieldname','inversion.cost_functions','size', [num_costfunc],'values',supportedcostfunctions())
                md = checkfield(md,'fieldname','inversion.cost_functions_coefficients','size',[md.mesh.numberofvertices, num_costfunc],'>=',0)
                md = checkfield(md,'fieldname','inversion.min_parameters','size',[md.mesh.numberofvertices, num_controls])
                md = checkfield(md,'fieldname','inversion.max_parameters','size',[md.mesh.numberofvertices, num_controls])
@@ -161,38 +162,38 @@
                        md = checkfield(md,'fieldname','inversion.vx_obs','size',[md.mesh.numberofvertices],'NaN',1,'Inf',1)
                        md = checkfield(md,'fieldname','inversion.vy_obs','size',[md.mesh.numberofvertices],'NaN',1,'Inf',1)
 
-               def marshall(self, md, fid):
+       def marshall(self,prefix,md,fid):
 
-                       yts=md.constants.yts;
-                       WriteData(fid,prefix,'object',self,'class','inversion','fieldname','iscontrol','format','Boolean')
-                       WriteData(fid,prefix,'name','md.inversion.type','data',1,'format','Integer')
-                       if not self.iscontrol:
-                               return
-                       WriteData(fid,prefix,'object',self,'class','inversion','fieldname','incomplete_adjoint','format','Boolean')
-                       WriteData(fid,prefix,'object',self,'class','inversion','fieldname','maxsteps','format','Integer')
-                       WriteData(fid,prefix,'object',self,'class','inversion','fieldname','maxiter','format','Integer')
-                       WriteData(fid,prefix,'object',self,'class','inversion','fieldname','fatol','format','Double')
-                       WriteData(fid,prefix,'object',self,'class','inversion','fieldname','frtol','format','Double')
-                       WriteData(fid,prefix,'object',self,'class','inversion','fieldname','gatol','format','Double')
-                       WriteData(fid,prefix,'object',self,'class','inversion','fieldname','grtol','format','Double')
-                       WriteData(fid,prefix,'object',self,'class','inversion','fieldname','gttol','format','Double')
-                       WriteData(fid,prefix,'object',self,'class','inversion','fieldname','algorithm','format','String')
-                       WriteData(fid,prefix,'object',self,'class','inversion','fieldname','cost_functions_coefficients','format','DoubleMat','mattype',1)
-                       WriteData(fid,prefix,'object',self,'class','inversion','fieldname','min_parameters','format','DoubleMat','mattype',3)
-                       WriteData(fid,prefix,'object',self,'class','inversion','fieldname','max_parameters','format','DoubleMat','mattype',3)
-                       WriteData(fid,prefix,'object',self,'class','inversion','fieldname','vx_obs','format','DoubleMat','mattype',1,'scale',1./yts)
-                       WriteData(fid,prefix,'object',self,'class','inversion','fieldname','vy_obs','format','DoubleMat','mattype',1,'scale',1./yts)
-                       WriteData(fid,prefix,'object',self,'class','inversion','fieldname','vz_obs','format','DoubleMat','mattype',1,'scale',1./yts)
-                       WriteData(fid,prefix,'object',self,'class','inversion','fieldname','thickness_obs','format','DoubleMat','mattype',1)
-                       WriteData(fid,prefix,'object',self,'class','inversion','fieldname','surface_obs','format','DoubleMat','mattype',1)
+               yts=md.constants.yts;
+               WriteData(fid,prefix,'object',self,'class','inversion','fieldname','iscontrol','format','Boolean')
+               WriteData(fid,prefix,'name','md.inversion.type','data',1,'format','Integer')
+               if not self.iscontrol:
+                       return
+               WriteData(fid,prefix,'object',self,'class','inversion','fieldname','incomplete_adjoint','format','Boolean')
+               WriteData(fid,prefix,'object',self,'class','inversion','fieldname','maxsteps','format','Integer')
+               WriteData(fid,prefix,'object',self,'class','inversion','fieldname','maxiter','format','Integer')
+               WriteData(fid,prefix,'object',self,'class','inversion','fieldname','fatol','format','Double')
+               WriteData(fid,prefix,'object',self,'class','inversion','fieldname','frtol','format','Double')
+               WriteData(fid,prefix,'object',self,'class','inversion','fieldname','gatol','format','Double')
+               WriteData(fid,prefix,'object',self,'class','inversion','fieldname','grtol','format','Double')
+               WriteData(fid,prefix,'object',self,'class','inversion','fieldname','gttol','format','Double')
+               WriteData(fid,prefix,'object',self,'class','inversion','fieldname','algorithm','format','String')
+               WriteData(fid,prefix,'object',self,'class','inversion','fieldname','cost_functions_coefficients','format','DoubleMat','mattype',1)
+               WriteData(fid,prefix,'object',self,'class','inversion','fieldname','min_parameters','format','DoubleMat','mattype',3)
+               WriteData(fid,prefix,'object',self,'class','inversion','fieldname','max_parameters','format','DoubleMat','mattype',3)
+               WriteData(fid,prefix,'object',self,'class','inversion','fieldname','vx_obs','format','DoubleMat','mattype',1,'scale',1./yts)
+               WriteData(fid,prefix,'object',self,'class','inversion','fieldname','vy_obs','format','DoubleMat','mattype',1,'scale',1./yts)
+               WriteData(fid,prefix,'object',self,'class','inversion','fieldname','vz_obs','format','DoubleMat','mattype',1,'scale',1./yts)
+               WriteData(fid,prefix,'object',self,'class','inversion','fieldname','thickness_obs','format','DoubleMat','mattype',1)
+               WriteData(fid,prefix,'object',self,'class','inversion','fieldname','surface_obs','format','DoubleMat','mattype',1)
 
-                       #process control parameters
-                       num_control_parameters = np.numel(self.control_parameters)
-                       WriteData(fid,prefix,'object',self,'fieldname','control_parameters','format','StringArray')
-                       WriteData(fid,prefix,'data',num_control_parameters,'name','md.inversion.num_control_parameters','format','Integer')
+               #process control parameters
+               num_control_parameters = np.size(self.control_parameters)
+               WriteData(fid,prefix,'object',self,'fieldname','control_parameters','format','StringArray')
+               WriteData(fid,prefix,'data',num_control_parameters,'name','md.inversion.num_control_parameters','format','Integer')
 
-                       #process cost functions
-                       num_cost_functions = np.size(self.cost_functions,2)
-                       data= marshallcostfunctions(self.cost_functions)
-                       WriteData(fid,prefix,'data',data,'name','md.inversion.cost_functions','format','StringArray')
-                       WriteData(fid,prefix,'data',num_cost_functions,'name','md.inversion.num_cost_functions','format','Integer')
+               #process cost functions
+               num_cost_functions = np.size(self.cost_functions)
+               data= marshallcostfunctions(self.cost_functions)
+               WriteData(fid,prefix,'data',data,'name','md.inversion.cost_functions','format','StringArray')
+               WriteData(fid,prefix,'data',num_cost_functions,'name','md.inversion.num_cost_functions','format','Integer')
Index: ../trunk-jpl/src/m/classes/SMBgemb.m
===================================================================
--- ../trunk-jpl/src/m/classes/SMBgemb.m        (revision 22266)
+++ ../trunk-jpl/src/m/classes/SMBgemb.m        (revision 22267)
@@ -357,7 +357,7 @@
                        dt=min(dtime);
             
                        WriteData(fid,prefix,'data',dt,'name','md.smb.dt','format','Double','scale',yts);
-            
+
                        % Check if smb_dt goes evenly into transient core time step
                        if (mod(md.timestepping.time_step,dt) >= 1e-10)
                 error('smb_dt/dt = %f. The number of SMB time steps in one transient core time step has to be an an integer',md.timestepping.time_step/dt);
Index: ../trunk-jpl/src/m/classes/plumebasalforcings.py
===================================================================
--- ../trunk-jpl/src/m/classes/plumebasalforcings.py    (nonexistent)
+++ ../trunk-jpl/src/m/classes/plumebasalforcings.py    (revision 22267)
@@ -0,0 +1,126 @@
+from fielddisplay import fielddisplay
+from checkfield import checkfield
+from WriteData import WriteData
+from project3d import project3d
+
+class plumebasalforcings(object):
+       '''
+       PLUME BASAL FORCINGS class definition
+
+               Usage:
+                       plumebasalforcings=plumebasalforcings()
+       '''
+
+       def __init__(self): # {{{
+               floatingice_melting_rate  = float('NaN')
+               groundedice_melting_rate  = float('NaN')
+               mantleconductivity        = float('NaN')
+               nusselt                   = float('NaN')
+               dtbg                      = float('NaN')
+               plumeradius               = float('NaN')
+               topplumedepth             = float('NaN')
+               bottomplumedepth          = float('NaN')
+               plumex                    = float('NaN')
+               plumey                    = float('NaN')
+               crustthickness            = float('NaN')
+               uppercrustthickness       = float('NaN')
+               uppercrustheat            = float('NaN')
+               lowercrustheat            = float('NaN')
+
+               self.setdefaultparameters()
+       #}}}
+
+       def __repr__(self): # {{{
+               print '   mantle plume basal melt parameterization:'
+
+               string="%s\n%s"%(string,fielddisplay(self,'groundedice_melting_rate','basal melting rate (positive if melting) [m/yr]'))
+               string="%s\n%s"%(string,fielddisplay(self,'floatingice_melting_rate','basal melting rate (positive if melting) [m/yr]'))
+               string="%s\n%s"%(string,fielddisplay(self,'mantleconductivity','mantle heat conductivity [W/m^3]'))
+               string="%s\n%s"%(string,fielddisplay(self,'nusselt','nusselt number, ratio of mantle to plume [1]'))
+               string="%s\n%s"%(string,fielddisplay(self,'dtbg','background temperature gradient [degree/m]'))
+               string="%s\n%s"%(string,fielddisplay(self,'plumeradius','radius of the mantle plume [m]'))
+               string="%s\n%s"%(string,fielddisplay(self,'topplumedepth','depth of the mantle plume top below the crust [m]'))
+               string="%s\n%s"%(string,fielddisplay(self,'bottomplumedepth','depth of the mantle plume base below the crust [m]'))
+               string="%s\n%s"%(string,fielddisplay(self,'plumex','x coordinate of the center of the plume [m]'))
+               string="%s\n%s"%(string,fielddisplay(self,'plumey','y coordinate of the center of the plume [m]'))
+               string="%s\n%s"%(string,fielddisplay(self,'crustthickness','thickness of the crust [m]'))
+               string="%s\n%s"%(string,fielddisplay(self,'uppercrustthickness','thickness of the upper crust [m]'))
+               string="%s\n%s"%(string,fielddisplay(self,'uppercrustheat','volumic heat of the upper crust [w/m^3]'))
+               string="%s\n%s"%(string,fielddisplay(self,'lowercrustheat','volumic heat of the lowercrust [w/m^3]'))
+
+               return string
+       #}}}
+
+       def initialize(self,md): #{{{
+               if np.all(np.isnan(self.groundedice_melting_rate)):
+                       self.groundedice_melting_rate=np.zeros((md.mesh.numberofvertices,))
+                       print '      no basalforcings.groundedice_melting_rate specified: values set as zero'
+               if np.all(np.isnan(self.floatingice_melting_rate)):
+                       self.floatingice_melting_rate=np.zeros((md.mesh.numberofvertices,))
+                       print '      no basalforcings.floatingice_melting_rate specified: values set as zero'
+               return
+       #}}}
+
+       def extrude(self,md): # {{{
+               self.groundedice_melting_rate=project3d(md,'vector',self.groundedice_melting_rate,'type','node','layer',1);
+               self.floatingice_melting_rate=project3d(md,'vector',self.floatingice_melting_rate,'type','node','layer',1);
+               return self
+       #}}}
+
+       def setdefaultparameters(self): # {{{
+               #default values for melting parameterization
+               self.mantleconductivity     = 2.2
+               self.nusselt                = 300
+               self.dtbg                   = 11/1000.
+               self.plumeradius            = 100000
+               self.topplumedepth          = 10000
+               self.bottomplumedepth       = 1050000
+               self.crustthickness         = 30000
+               self.uppercrustthickness    = 14000
+               self.uppercrustheat         = 1.7*10**-6
+               self.lowercrustheat         = 0.4*10**-6
+               return self
+       #}}}
+
+       def checkconsistency(self,md,solution,analyses):    # {{{
+               if 'MasstransportAnalysis' in analyses and not (solution == 'TransientSolution' and md.transient.ismasstransport==0):
+                       md = checkfield(md,'fieldname','basalforcings.groundedice_melting_rate','NaN',1,'timeseries',1)
+                       md = checkfield(md,'fieldname','basalforcings.floatingice_melting_rate','NaN',1,'timeseries',1)
+               if 'BalancethicknessAnalysis' in analyses:
+                       md = checkfield(md,'fieldname','basalforcings.groundedice_melting_rate','NaN',1,'size',[md.mesh.numberofvertices])
+                       md = checkfield(md,'fieldname','basalforcings.floatingice_melting_rate','NaN',1,'size',[md.mesh.numberofvertices])
+               if 'ThermalAnalysis' in analyses and not (solution == 'TransientSolution' and md.transient.isthermal==0):
+                       md = checkfield(md,'fieldname','basalforcings.groundedice_melting_rate','NaN',1,'timeseries',1)
+                       md = checkfield(md,'fieldname','basalforcings.floatingice_melting_rate','NaN',1,'timeseries',1)
+                       md = checkfield(md,'fieldname','basalforcings.mantleconductivity','>=',0,'numel',1)
+                       md = checkfield(md,'fieldname','basalforcings.nusselt','>',0,'numel',1)
+                       md = checkfield(md,'fieldname','basalforcings.dtbg','>',0,'numel',1)
+                       md = checkfield(md,'fieldname','basalforcings.topplumedepth','>',0,'numel',1)
+                       md = checkfield(md,'fieldname','basalforcings.bottomplumedepth','>',0,'numel',1)
+                       md = checkfield(md,'fieldname','basalforcings.plumex','numel',1)
+                       md = checkfield(md,'fieldname','basalforcings.plumey','numel',1)
+                       md = checkfield(md,'fieldname','basalforcings.crustthickness','>',0,'numel',1)
+                       md = checkfield(md,'fieldname','basalforcings.uppercrustthickness','>',0,'numel',1)
+                       md = checkfield(md,'fieldname','basalforcings.uppercrustheat','>',0,'numel',1)
+                       md = checkfield(md,'fieldname','basalforcings.lowercrustheat','>',0,'numel',1)
+               return md
+       # }}}
+       def marshall(self,prefix,md,fid):    # {{{
+               yts=md.constants.yts
+
+               WriteData(fid,prefix,'name','md.basalforcings.model','data',4,'format','Integer')
+               WriteData(fid,prefix,'object',self,'fieldname','floatingice_melting_rate','format','DoubleMat','name','md.basalforcings.floatingice_melting_rate','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts)
+               WriteData(fid,prefix,'object',self,'fieldname','groundedice_melting_rate','format','DoubleMat','name','md.basalforcings.groundedice_melting_rate','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts)
+               WriteData(fid,prefix,'object',self,'fieldname','mantleconductivity','format','Double')
+               WriteData(fid,prefix,'object',self,'fieldname','nusselt','format','Double')
+               WriteData(fid,prefix,'object',self,'fieldname','dtbg','format','Double')
+               WriteData(fid,prefix,'object',self,'fieldname','plumeradius','format','Double')
+               WriteData(fid,prefix,'object',self,'fieldname','topplumedepth','format','Double')
+               WriteData(fid,prefix,'object',self,'fieldname','bottomplumedepth','format','Double')
+               WriteData(fid,prefix,'object',self,'fieldname','plumex','format','Double')
+               WriteData(fid,prefix,'object',self,'fieldname','plumey','format','Double')
+               WriteData(fid,prefix,'object',self,'fieldname','crustthickness','format','Double')
+               WriteData(fid,prefix,'object',self,'fieldname','uppercrustthickness','format','Double')
+               WriteData(fid,prefix,'object',self,'fieldname','uppercrustheat','format','Double')
+               WriteData(fid,prefix,'object',self,'fieldname','lowercrustheat','format','Double')
+       # }}}
Index: ../trunk-jpl/test/NightlyRun/test701.py
===================================================================
--- ../trunk-jpl/test/NightlyRun/test701.py     (nonexistent)
+++ ../trunk-jpl/test/NightlyRun/test701.py     (revision 22267)
@@ -0,0 +1,75 @@
+#Test Name: FlowbandFSshelf
+import numpy as np
+from scipy.interpolate import interp1d
+from model import *
+from solve import *
+from setflowequation import *
+from bamgflowband import *
+from paterson import *
+
+x = np.arange(1,3001,100).T
+h = np.linspace(1000,300,np.size(x)).T
+b = -917./1023.*h
+
+md = bamgflowband(model(),x,b+h,b,'hmax',80.)
+print md.mesh.numberofvertices
+
+#Geometry          #interp1d returns a function to be called on md.mesh.x
+md.geometry.surface = interp1d(x,b+h)(md.mesh.x)
+md.geometry.base = interp1d(x,b)(md.mesh.x)
+md.geometry.thickness = md.geometry.surface-md.geometry.base
+
+#mask
+md.mask.ice_levelset = -np.ones((md.mesh.numberofvertices,))
+md.mask.ice_levelset[np.where(md.mesh.vertexflags(2))] = 0.
+md.mask.groundedice_levelset = np.zeros((md.mesh.numberofvertices,)) - 0.5
+
+#materials
+md.initialization.temperature = (273.-20.)*np.ones((md.mesh.numberofvertices,))
+md.materials.rheology_B = paterson(md.initialization.temperature)
+md.materials.rheology_n = 3.*np.ones((md.mesh.numberofelements,))
+
+#friction
+md.friction.coefficient = np.zeros((md.mesh.numberofvertices,))
+md.friction.coefficient[np.where(md.mesh.vertexflags(1))] = 20.
+md.friction.p = np.ones((md.mesh.numberofelements,))
+md.friction.q = np.ones((md.mesh.numberofelements,))
+
+#Boundary conditions
+md.stressbalance.referential = float('NaN')*np.ones((md.mesh.numberofvertices,6))
+md.stressbalance.loadingforce = 0. * np.ones((md.mesh.numberofvertices,3))
+md.stressbalance.spcvx = float('NaN') * np.ones((md.mesh.numberofvertices,))
+md.stressbalance.spcvy = float('NaN') * np.ones((md.mesh.numberofvertices,))
+md.stressbalance.spcvz = float('NaN') * np.ones((md.mesh.numberofvertices,))
+md.stressbalance.spcvx[np.where(md.mesh.vertexflags(4))] = 0.
+md.stressbalance.spcvy[np.where(md.mesh.vertexflags(4))] = 0.
+
+#Misc
+md = setflowequation(md,'FS','all')
+md.stressbalance.abstol = float('NaN')
+#md.stressbalance.reltol = 10**-16
+md.stressbalance.FSreconditioning = 1.
+md.stressbalance.maxiter = 20
+md.flowequation.augmented_lagrangian_r = 10000.
+md.miscellaneous.name  =  'test701'
+md.verbose = verbose('convergence',True)
+md.cluster = generic('np',2)
+
+#Go solve
+field_names = []
+field_tolerances = []
+field_values = []
+#md.initialization.pressure = md.constants.g*md.materials.rho_ice*(md.geometry.surface-md.mesh.y)
+for i in ['MINI','MINIcondensed','TaylorHood','LATaylorHood','CrouzeixRaviart','LACrouzeixRaviart']:
+       print ' '
+       print '======Testing ' +i+ ' Full-Stokes Finite element====='
+       md.flowequation.fe_FS = i
+       md = solve(md,'Stressbalance')
+       field_names = field_names + [['Vx'+i],['Vy'+i],['Vel'+i],['Pressure'+i]]
+       field_tolerances = field_tolerances + [9e-5,9e-5,9e-5,1e-10]
+       field_values = field_values + [
+               md.results.StressbalanceSolution.Vx,
+               md.results.StressbalanceSolution.Vy,
+               md.results.StressbalanceSolution.Vel,
+               md.results.StressbalanceSolution.Pressure
+               ]
Index: ../trunk-jpl/test/NightlyRun/test703.py
===================================================================
--- ../trunk-jpl/test/NightlyRun/test703.py     (nonexistent)
+++ ../trunk-jpl/test/NightlyRun/test703.py     (revision 22267)
@@ -0,0 +1,160 @@
+#Test Name: FlowbandFSsheetshelfTrans
+import numpy as np
+from scipy.interpolate import interp1d
+import copy
+from model import  * 
+from socket import gethostname
+from triangle import  * 
+from meshconvert import  * 
+from setmask import  * 
+from parameterize import  * 
+from setflowequation import  * 
+from solve import  *
+from NowickiProfile import *
+from bamgflowband import *
+from paterson import *
+
+#mesh parameters
+x  = np.arange(-5,5.5,.5).T
+[b,h,sea] = NowickiProfile(x)
+x = x * 10**3
+h = h * 10**3
+b = (b-sea) * 10**3
+
+#mesh domain
+md = bamgflowband(model(),x,b+h,b,'hmax',150)
+print md.mesh.numberofvertices
+
+#geometry
+md.geometry.surface = interp1d(x,b+h)(md.mesh.x)
+md.geometry.base = interp1d(x,b)(md.mesh.x)
+md.geometry.thickness = md.geometry.surface - md.geometry.base
+
+#mask
+md.mask.ice_levelset = -np.ones((md.mesh.numberofvertices,))
+md.mask.ice_levelset[np.where(md.mesh.vertexflags(2))] = 0
+md.mask.groundedice_levelset = np.zeros((md.mesh.numberofvertices,)) - 0.5
+
+#materials
+md.initialization.temperature = (273. - 20.) * np.ones((md.mesh.numberofvertices,))
+md.materials.rheology_B = paterson(md.initialization.temperature)
+md.materials.rheology_n = 3 * np.ones((md.mesh.numberofelements,))
+
+#friction
+md.friction.coefficient = np.zeros((md.mesh.numberofvertices,))
+md.friction.coefficient[np.where(md.mesh.vertexflags(1))] = 20
+md.friction.p = np.ones((md.mesh.numberofelements,))
+md.friction.q = np.ones((md.mesh.numberofelements,))
+
+#boundary conditions
+md.stressbalance.spcvx = float('NaN') * np.ones((md.mesh.numberofvertices,))
+md.stressbalance.spcvy = float('NaN') * np.ones((md.mesh.numberofvertices,))
+md.stressbalance.spcvz = float('NaN') * np.ones((md.mesh.numberofvertices,))
+md.stressbalance.referential = float('NaN') * np.ones((md.mesh.numberofvertices,6))
+md.stressbalance.loadingforce = 0 * np.ones((md.mesh.numberofvertices,3))
+md.stressbalance.spcvx[np.where(md.mesh.vertexflags(4))] = 800.
+md.stressbalance.spcvy[np.where(md.mesh.vertexflags(4))] = 0.
+
+#print md.stressbalance.referential
+#print md.stressbalance.loadingforce
+#print md.stressbalance.spcvx
+#print md.stressbalance.spcvy
+#print md.stressbalance.spcvz
+#print np.shape(md.stressbalance.referential)
+#print np.shape(md.stressbalance.loadingforce)
+#print np.shape(md.stressbalance.spcvx)
+#print np.shape(md.stressbalance.spcvy)
+#print np.shape(md.stressbalance.spcvz)
+
+#Misc
+md = setflowequation(md,'FS','all')
+md.flowequation.fe_FS = 'TaylorHood'
+md.stressbalance.abstol = float('NaN')
+md.miscellaneous.name = 'test703'
+
+#Transient settings
+md.timestepping.time_step = 0.000001
+md.timestepping.final_time = 0.000005
+md.stressbalance.shelf_dampening = 1.
+md.smb.mass_balance = np.zeros((md.mesh.numberofvertices,))
+md.basalforcings.groundedice_melting_rate = np.zeros((md.mesh.numberofvertices,))
+md.basalforcings.floatingice_melting_rate = np.zeros((md.mesh.numberofvertices,))
+md.basalforcings.geothermalflux = np.zeros((md.mesh.numberofvertices,))
+posb = np.intersect1d(np.where(md.mesh.x > 0.), np.where(md.mesh.vertexonbase))
+md.basalforcings.groundedice_melting_rate[posb] = 18.
+md.basalforcings.floatingice_melting_rate[posb] = 18.
+md.initialization.vx = np.zeros((md.mesh.numberofvertices,))
+md.initialization.vy = np.zeros((md.mesh.numberofvertices,))
+md.initialization.pressure = np.zeros((md.mesh.numberofvertices,))
+md.masstransport.spcthickness = float('NaN') * np.ones((md.mesh.numberofvertices,))
+md.thermal.spctemperature = float('NaN') * np.ones((md.mesh.numberofvertices,))
+md.transient.isthermal = 0
+md.masstransport.isfreesurface = 1
+
+#Go solve
+md.cluster = generic('np',3)
+md.stressbalance.shelf_dampening = 1
+md1 = solve(md,'Transient')
+
+md.stressbalance.shelf_dampening = 0
+md = solve(md,'Transient')
+
+#Fields and tolerances to track changes
+field_names      = [
+       'Vx1','Vy1','Vel1','Pressure1','Bed1','Surface1','Thickness1',
+       'Vx2','Vy2','Vel2','Pressure2','Bed2','Surface2','Thickness2',
+       'Vx3','Vy3','Vel3','Pressure3','Bed3','Surface3','Thickness3',
+       'Vx1_damp','Vy1_damp','Vel1_damp','Pressure1_damp','Bed1_damp','Surface1_damp','Thickness1_damp',
+       'Vx2_damp','Vy2_damp','Vel2_damp','Pressure2_damp','Bed2_damp','Surface2_damp','Thickness2_damp',
+       'Vx3_damp','Vy3_damp','Vel3_damp','Pressure3_damp','Bed3_damp','Surface3_damp','Thickness3_damp']
+field_tolerances = [
+       2e-08,2e-08,2e-08,1e-08,1e-10,1e-10,1e-10,
+       2e-08,2e-08,2e-08,1e-08,1e-10,1e-10,1e-10,
+       2e-08,2e-08,2e-08,1e-08,1e-10,1e-10,1e-10,
+       5e-08,5e-08,5e-08,1e-08,1e-10,1e-10,1e-10,
+       5e-08,5e-08,5e-08,1e-08,1e-10,1e-10,1e-10,
+       5e-08,5e-08,5e-08,1e-08,1e-10,1e-10,1e-10]
+field_values = [
+       md.results.TransientSolution[0].Vx,
+       md.results.TransientSolution[0].Vy,
+       md.results.TransientSolution[0].Vel,
+       md.results.TransientSolution[0].Pressure,
+       md.results.TransientSolution[0].Base,
+       md.results.TransientSolution[0].Surface,
+       md.results.TransientSolution[0].Thickness,
+       md.results.TransientSolution[1].Vx,
+       md.results.TransientSolution[1].Vy,
+       md.results.TransientSolution[1].Vel,
+       md.results.TransientSolution[1].Pressure,
+       md.results.TransientSolution[1].Base,
+       md.results.TransientSolution[1].Surface,
+       md.results.TransientSolution[1].Thickness,
+       md.results.TransientSolution[2].Vx,
+       md.results.TransientSolution[2].Vy,
+       md.results.TransientSolution[2].Vel,
+       md.results.TransientSolution[2].Pressure,
+       md.results.TransientSolution[2].Base,
+       md.results.TransientSolution[2].Surface,
+       md.results.TransientSolution[2].Thickness,
+       md1.results.TransientSolution[0].Vx,
+       md1.results.TransientSolution[0].Vy,
+       md1.results.TransientSolution[0].Vel,
+       md1.results.TransientSolution[0].Pressure,
+       md1.results.TransientSolution[0].Base,
+       md1.results.TransientSolution[0].Surface,
+       md1.results.TransientSolution[0].Thickness,
+       md1.results.TransientSolution[1].Vx,
+       md1.results.TransientSolution[1].Vy,
+       md1.results.TransientSolution[1].Vel,
+       md1.results.TransientSolution[1].Pressure,
+       md1.results.TransientSolution[1].Base,
+       md1.results.TransientSolution[1].Surface,
+       md1.results.TransientSolution[1].Thickness,
+       md1.results.TransientSolution[2].Vx,
+       md1.results.TransientSolution[2].Vy,
+       md1.results.TransientSolution[2].Vel,
+       md1.results.TransientSolution[2].Pressure,
+       md1.results.TransientSolution[2].Base,
+       md1.results.TransientSolution[2].Surface,
+       md1.results.TransientSolution[2].Thickness,
+       ]
Index: ../trunk-jpl/test/NightlyRun/test293.py
===================================================================
--- ../trunk-jpl/test/NightlyRun/test293.py     (nonexistent)
+++ ../trunk-jpl/test/NightlyRun/test293.py     (revision 22267)
@@ -0,0 +1,60 @@
+#Test Name: SquareShelfTranSSA2dMismipFloatingMeltParam
+import numpy as np
+from model import *
+from socket import gethostname
+from triangle import *
+from setmask import *
+from parameterize import *
+from setflowequation import *
+from solve import *
+from mismipbasalforcings import *
+
+md = triangle(model(),'../Exp/Square.exp',150000.)
+md = setmask(md,'all','')
+md = parameterize(md,'../Par/SquareShelf.py')
+md = setflowequation(md,'SSA','all')
+md.cluster = generic('name',gethostname(),'np',3)
+md.constants.yts = 365.2422 * 24. * 3600.
+md.basalforcings = mismipbasalforcings()
+md.basalforcings = md.basalforcings.initialize(md)
+md.transient.isgroundingline = 1
+md.geometry.bed = min(md.geometry.base) * np.ones(md.mesh.numberofvertices,)
+md.transient.requested_outputs = ['default','BasalforcingsFloatingiceMeltingRate']
+print md.constants.yts
+
+md = solve(md,'Transient')
+
+#Fields and tolerances to track changes
+field_names     =['Vx1','Vy1','Vel1','Pressure1','Bed1','Surface1','Thickness1','BasalforcingsFloatingiceMeltingRate1',
+       'Vx2','Vy2','Vel2','Pressure2','Bed2','Surface2','Thickness2','BasalforcingsFloatingiceMeltingRate2',
+       'Vx3','Vy3','Vel3','Pressure3','Bed3','Surface3','Thickness3','BasalforcingsFloatingiceMeltingRate3']
+field_tolerances = [
+       1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,2e-13,
+       1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,2e-13,
+       1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,2e-13]
+field_values = [
+       md.results.TransientSolution[0].Vx,
+       md.results.TransientSolution[0].Vy,
+       md.results.TransientSolution[0].Vel,
+       md.results.TransientSolution[0].Pressure,
+       md.results.TransientSolution[0].Base,
+       md.results.TransientSolution[0].Surface,
+       md.results.TransientSolution[0].Thickness,
+       md.results.TransientSolution[0].BasalforcingsFloatingiceMeltingRate,
+       md.results.TransientSolution[1].Vx,
+       md.results.TransientSolution[1].Vy,
+       md.results.TransientSolution[1].Vel,
+       md.results.TransientSolution[1].Pressure,
+       md.results.TransientSolution[1].Base,
+       md.results.TransientSolution[1].Surface,
+       md.results.TransientSolution[1].Thickness,
+       md.results.TransientSolution[1].BasalforcingsFloatingiceMeltingRate,
+       md.results.TransientSolution[2].Vx,
+       md.results.TransientSolution[2].Vy,
+       md.results.TransientSolution[2].Vel,
+       md.results.TransientSolution[2].Pressure,
+       md.results.TransientSolution[2].Base,
+       md.results.TransientSolution[2].Surface,
+       md.results.TransientSolution[2].Thickness,
+       md.results.TransientSolution[2].BasalforcingsFloatingiceMeltingRate,
+       ]
Index: ../trunk-jpl/src/m/mesh/bamgflowband.py
===================================================================
--- ../trunk-jpl/src/m/mesh/bamgflowband.py     (nonexistent)
+++ ../trunk-jpl/src/m/mesh/bamgflowband.py     (revision 22267)
@@ -0,0 +1,47 @@
+import numpy as  np
+from model import *
+from collections import OrderedDict
+from bamg import *
+from mesh2dvertical import *
+
+def bamgflowband(md,x,surf,base,*args):
+       """
+       BAMGFLOWBAND - create flowband mesh with bamg
+
+       Usage:
+               md=bamgflowband(md,x,surf,base,OPTIONS)
+
+               surf and bed are the surface elevation and base for each x provided
+               x must be increasing
+               OPTIONS are bamg options
+
+       Example:
+               x =np.arrange(1,3001,100)
+               h=linspace(1000,300,numel(x))
+               b=-917/1023*h
+               md=bamgflowband(model,b+h,b,'hmax',80,'vertical',1,'Markers',m)
+       """
+
+       #Write expfile with domain outline
+       A = OrderedDict()
+       A.x = np.concatenate((x,np.flipud(x),[x[0]]))
+       A.y = np.concatenate((base,np.flipud(surf),[base[0]]))
+       A.nods = np.size(A.x)
+
+       #markers:
+       m                               = np.ones((np.size(A.x)-1,))    # base        = 1
+       m[np.size(x) - 1]                       = 2                     # right side  = 2
+       m[np.size(x):2 * np.size(x) - 1]        = 3                     # top surface = 3
+       m[2 * np.size(x) - 1]                   = 4                     # left side   = 4
+
+       #mesh domain
+       md = bamg(model(),'domain',[A],'vertical',1,'Markers',m,*args)
+       #print md.mesh.numberofvertices
+
+       #Deal with vertices on bed
+       md.mesh.vertexonbase = np.zeros((md.mesh.numberofvertices,))
+       md.mesh.vertexonbase[np.where(md.mesh.vertexflags(1))] = 1
+       md.mesh.vertexonsurface = np.zeros((md.mesh.numberofvertices,))
+       md.mesh.vertexonsurface[np.where(md.mesh.vertexflags(3))] = 1
+
+       return md
Index: ../trunk-jpl/src/m/mesh/bamg.py
===================================================================
--- ../trunk-jpl/src/m/mesh/bamg.py     (revision 22266)
+++ ../trunk-jpl/src/m/mesh/bamg.py     (revision 22267)
@@ -1,6 +1,6 @@
 import os.path
 import numpy as  np
-from mesh2d import mesh2d
+from mesh2dvertical import *
 from collections import OrderedDict
 from pairoptions import pairoptions
 from bamggeom import bamggeom
@@ -79,9 +79,12 @@
 
                #Check that file exists
                domainfile=options.getfieldvalue('domain')
-               if not os.path.exists(domainfile):
-                       raise IOError("bamg error message: file '%s' not found" % domainfile)
-               domain=expread(domainfile)
+               if type(domainfile) == str:
+                       if not os.path.exists(domainfile):
+                               raise IOError("bamg error message: file '%s' not found" % domainfile)
+                       domain=expread(domainfile)
+               else:
+                       domain=domainfile
 
                #Build geometry 
                count=0
@@ -88,18 +91,18 @@
                for i,domaini in enumerate(domain):
 
                        #Check that the domain is closed
-                       if (domaini['x'][0]!=domaini['x'][-1] or domaini['y'][0]!=domaini['y'][-1]):
+                       if (domaini.x[0]!=domaini.x[-1] or domaini.y[0]!=domaini.y[-1]):
                                raise RuntimeError("bamg error message: all contours provided in ''domain'' should be closed")
 
                        #Checks that all holes are INSIDE the principle domain outline
                        if i:
-                               flags=ContourToNodes(domaini['x'],domaini['y'],domainfile,0)[0]
+                               flags=ContourToNodes(domaini.x,domaini.y,domainfile,0)[0]
                                if np.any(np.logical_not(flags)):
                                        raise RuntimeError("bamg error message: All holes should be strictly inside the principal domain")
 
                        #Add all points to bamg_geometry
-                       nods=domaini['nods']-1    #the domain are closed 0=end
-                       bamg_geometry.Vertices=np.vstack((bamg_geometry.Vertices,np.vstack((domaini['x'][0:nods],domaini['y'][0:nods],np.ones((nods)))).T))
+                       nods=domaini.nods-1    #the domain are closed 0=end
+                       bamg_geometry.Vertices=np.vstack((bamg_geometry.Vertices,np.vstack((domaini.x[0:nods],domaini.y[0:nods],np.ones((nods)))).T))
                        bamg_geometry.Edges   =np.vstack((bamg_geometry.Edges,np.vstack((np.arange(count+1,count+nods+1),np.hstack((np.arange(count+2,count+nods+1),count+1)),1.*np.ones((nods)))).T))
                         if i:
                             bamg_geometry.SubDomains=np.vstack((bamg_geometry.SubDomains,[2,count+1,1,-subdomain_ref]))
@@ -115,6 +118,12 @@
                        #update counter
                        count+=nods
 
+               if options.getfieldvalue('vertical',0):
+                       if np.size(options.getfieldvalue('Markers',[]))!=np.size(bamg_geometry.Edges,0):
+                               raise RuntimeError('for 2d vertical mesh, ''Markers'' option is required, and should be of size ' + str(np.size(bamg_geometry.Edges,0)))
+                       if np.size(options.getfieldvalue('Markers',[]))==np.size(bamg_geometry.Edges,0):
+                               bamg_geometry.Edges[:,2]=options.getfieldvalue('Markers')
+
                #take care of rifts
                if options.exist('rifts'):
 
@@ -322,23 +331,64 @@
        bamgmesh_out,bamggeom_out=BamgMesher(bamg_mesh.__dict__,bamg_geometry.__dict__,bamg_options)
 
        # plug results onto model
+       if options.getfieldvalue('vertical',0):
+               md.mesh=mesh2dvertical()
+               md.mesh.x=bamgmesh_out['Vertices'][:,0].copy()
+               md.mesh.y=bamgmesh_out['Vertices'][:,1].copy()
+               md.mesh.elements=bamgmesh_out['Triangles'][:,0:3].astype(int)
+               md.mesh.edges=bamgmesh_out['IssmEdges'].astype(int)
+               md.mesh.segments=bamgmesh_out['IssmSegments'][:,0:3].astype(int)
+               md.mesh.segmentmarkers=bamgmesh_out['IssmSegments'][:,3].astype(int)
+
+               #Fill in rest of fields:
+               md.mesh.numberofelements=np.size(md.mesh.elements,axis=0)
+               md.mesh.numberofvertices=np.size(md.mesh.x)
+               md.mesh.numberofedges=np.size(md.mesh.edges,axis=0)
+               md.mesh.vertexonboundary=np.zeros(md.mesh.numberofvertices,bool)
+               md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1]=True
+
+               #Now, build the connectivity tables for this mesh. Doubled in matlab for some reason
+               md.mesh.vertexonboundary=np.zeros(md.mesh.numberofvertices,)
+               md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1]=1
+
+       elif options.getfieldvalue('3dsurface',0):
+               md.mesh=mesh3dsurface()
+               md.mesh.x=bamgmesh_out['Vertices'][:,0].copy()
+               md.mesh.y=bamgmesh_out['Vertices'][:,1].copy()
+               md.mesh.z=md.mesh.x
+               md.mesh.z[:]=0
+               md.mesh.elements=bamgmesh_out['Triangles'][:,0:3].astype(int)
+               md.mesh.edges=bamgmesh_out['IssmEdges'].astype(int)
+               md.mesh.segments=bamgmesh_out['IssmSegments'][:,0:3].astype(int)
+               md.mesh.segmentmarkers=bamgmesh_out['IssmSegments'][:,3].astype(int)
+
+               #Fill in rest of fields:
+               md.mesh.numberofelements=np.size(md.mesh.elements,axis=0)
+               md.mesh.numberofvertices=np.size(md.mesh.x)
+               md.mesh.numberofedges=np.size(md.mesh.edges,axis=0)
+               md.mesh.vertexonboundary=np.zeros(md.mesh.numberofvertices,bool)
+               md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1]=True
+
+       else:
+               md.mesh=mesh2d()
+               md.mesh.x=bamgmesh_out['Vertices'][:,0].copy()
+               md.mesh.y=bamgmesh_out['Vertices'][:,1].copy()
+               md.mesh.elements=bamgmesh_out['Triangles'][:,0:3].astype(int)
+               md.mesh.edges=bamgmesh_out['IssmEdges'].astype(int)
+               md.mesh.segments=bamgmesh_out['IssmSegments'][:,0:3].astype(int)
+               md.mesh.segmentmarkers=bamgmesh_out['IssmSegments'][:,3].astype(int)
+
+               #Fill in rest of fields:
+               md.mesh.numberofelements=np.size(md.mesh.elements,axis=0)
+               md.mesh.numberofvertices=np.size(md.mesh.x)
+               md.mesh.numberofedges=np.size(md.mesh.edges,axis=0)
+               md.mesh.vertexonboundary=np.zeros(md.mesh.numberofvertices,bool)
+               md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1]=True
+
+       #Bamg private fields
        md.private.bamg=OrderedDict()
        md.private.bamg['mesh']=bamgmesh(bamgmesh_out)
        md.private.bamg['geometry']=bamggeom(bamggeom_out)
-       md.mesh = mesh2d()
-       md.mesh.x=bamgmesh_out['Vertices'][:,0].copy()
-       md.mesh.y=bamgmesh_out['Vertices'][:,1].copy()
-       md.mesh.elements=bamgmesh_out['Triangles'][:,0:3].astype(int)
-       md.mesh.edges=bamgmesh_out['IssmEdges'].astype(int)
-       md.mesh.segments=bamgmesh_out['IssmSegments'][:,0:3].astype(int)
-       md.mesh.segmentmarkers=bamgmesh_out['IssmSegments'][:,3].astype(int)
-
-       #Fill in rest of fields:
-       md.mesh.numberofelements=np.size(md.mesh.elements,axis=0)
-       md.mesh.numberofvertices=np.size(md.mesh.x)
-       md.mesh.numberofedges=np.size(md.mesh.edges,axis=0)
-       md.mesh.vertexonboundary=np.zeros(md.mesh.numberofvertices,bool)
-       md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1]=True
        md.mesh.elementconnectivity=md.private.bamg['mesh'].ElementConnectivity
        md.mesh.elementconnectivity[np.nonzero(np.isnan(md.mesh.elementconnectivity))]=0
        md.mesh.elementconnectivity=md.mesh.elementconnectivity.astype(int)