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

../trunk-jpl/test/Par/SquareThermal.py

@@ -25 +25 @@
 
 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 +35 @@
 
 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 +45 @@
 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

../trunk-jpl/test/Par/ISMIPE.par

@@ -13 +13 @@
         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;

../trunk-jpl/test/Par/ISMIPE.py

@@ -12 +12 @@
         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

../trunk-jpl/test/NightlyRun/test2425.py

@@ +1 @@
+#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]
+

../trunk-jpl/test/NightlyRun/test342.py

@@ +1 @@
+#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,
+        ]

../trunk-jpl/test/NightlyRun/test261.py

@@ +1 @@
+#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,
+        ]

../trunk-jpl/test/NightlyRun/test441.py

@@ +1 @@
+#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,
+        ]

../trunk-jpl/test/NightlyRun/test343.py

@@ +1 @@
+#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
+        ]

../trunk-jpl/test/NightlyRun/test540.py

@@ +1 @@
+#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,
+        ]

../trunk-jpl/test/NightlyRun/test442.py

@@ +1 @@
+#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,
+        ]

../trunk-jpl/test/NightlyRun/test344.py

@@ +1 @@
+#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,
+        ]

../trunk-jpl/test/NightlyRun/test460.py

@@ +1 @@
+#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,
+                        ]

../trunk-jpl/test/NightlyRun/test461.py

@@ +1 @@
+#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,
+           ]

../trunk-jpl/test/NightlyRun/test435.py

@@ +1 @@
+#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,
+        ]

../trunk-jpl/test/NightlyRun/test462.py

@@ +1 @@
+#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,
+        ]

../trunk-jpl/test/NightlyRun/test436.py

@@ +1 @@
+#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,
+                ]
+

../trunk-jpl/test/NightlyRun/test463.py

@@ +1 @@
+#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,
+        ]

../trunk-jpl/test/NightlyRun/test437.py

@@ +1 @@
+#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,
+        ]

../trunk-jpl/test/NightlyRun/test464.py

@@ +1 @@
+#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,
+        ]

../trunk-jpl/test/NightlyRun/test438.py

@@ +1 @@
+#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
+        ]

../trunk-jpl/test/NightlyRun/test465.py

@@ +1 @@
+#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,
+        ]

../trunk-jpl/test/NightlyRun/test439.py

@@ +1 @@
+#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
+        ]

../trunk-jpl/test/NightlyRun/test340.py

@@ +1 @@
+#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
+]

../trunk-jpl/test/NightlyRun/test430.py

@@ +1 @@
+#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,
+        ]

../trunk-jpl/test/NightlyRun/test808.py

@@ +1 @@
+#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
+        ]

../trunk-jpl/test/NightlyRun/test2424.py

@@ +1 @@
+#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))
+

../trunk-jpl/test/NightlyRun/test260.py

@@ +1 @@
+#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,
+        ]

../trunk-jpl/test/NightlyRun/test350.py

@@ +1 @@
+#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
+        ]
+

../trunk-jpl/test/NightlyRun/test243.py

@@ +1 @@
+#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],
+        ]

../trunk-jpl/src/m/solve/WriteData.py

@@ -50 +50 @@
         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)))

../trunk-jpl/src/m/solve/parseresultsfromdisk.py

@@ -174 +174 @@
                 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 +188 @@
                         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

../trunk-jpl/src/m/geometry/NowickiProfile.py

@@ +1 @@
+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]

../trunk-jpl/src/m/consistency/checkfield.py

@@ -83 +83 @@
         #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 +100 @@
                         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 +107 @@
                         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 +130 @@
 
         #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 +158 @@
         #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 +194 @@
 
         #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 +228 @@
 
         return md
 
+

../trunk-jpl/src/m/mech/newforcing.py

@@ +1 @@
+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

../trunk-jpl/src/m/dev/devpath.py

@@ -30 +30 @@
 
 #Manual imports for commonly used functions
 from plotmodel import plotmodel
+from runme import runme
 
 #c = get_ipython().config
 #c.InteractiveShellApp.exec_lines = []

../trunk-jpl/src/m/print/printmodel.py

@@ +1 @@
+
+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

../trunk-jpl/src/m/classes/thermal.py

@@ -100 +100 @@
                 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):

../trunk-jpl/src/m/classes/calvingminthickness.py

@@ +1 @@
+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)
+        # }}}

../trunk-jpl/src/m/classes/calvingdev.py

@@ +1 @@
+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)
+        # }}}

../trunk-jpl/src/m/classes/mismipbasalforcings.py

@@ -39 +39 @@
     #}}}
     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 +86 @@
             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)

../trunk-jpl/src/m/classes/SMBgradientsela.py

@@ -8 +8 @@
         SMBgradientsela Class definition
 
            Usage:
-              SMBgradientsela=SMBgradientsela();
+              SMBgradientsela=SMBgradientsela()
         """
 
         def __init__(self): # {{{
@@ -18 +18 @@
                 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 +47 @@
 
                 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 +61 @@
                         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):    # {{{

../trunk-jpl/src/m/classes/constants.py

@@ -11 +11 @@
         """
 
         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 +48 @@
         #}}}
         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
         # }}}

../trunk-jpl/src/m/classes/hydrologysommers.py

@@ +1 @@
+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')
+        # }}}

../trunk-jpl/src/m/classes/flowequation.py

@@ -91 +91 @@
                 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 +103 @@
                 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))

../trunk-jpl/src/m/classes/basalforcings.py

@@ -44 +44 @@
                 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
         #}}}

../trunk-jpl/src/m/classes/calvingvonmises.py

@@ +1 @@
+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)
+        # }}}

../trunk-jpl/src/m/classes/transient.py

@@ -11 +11 @@
         """
 
         def __init__(self): # {{{
-                self.issmb   = False
+                self.issmb             = False
                 self.ismasstransport   = False
                 self.isstressbalance   = False
                 self.isthermal         = False
@@ -22 +22 @@
                 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 +80 @@
                 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

../trunk-jpl/src/m/classes/mesh2dvertical.py

@@ +1 @@
+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')
+        # }}}

../trunk-jpl/src/m/classes/matenhancedice.py

@@ +1 @@
+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')
+        # }}}

../trunk-jpl/src/m/classes/amr.py

@@ -27 +27 @@
         self.deviatoricerror_resolution= 0.
         self.deviatoricerror_threshold = 0.
         self.deviatoricerror_maximum    = 0.
+        self.restart=0.
         #set defaults
         self.setdefaultparameters()
     #}}}
@@ -47 +48 @@
         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 +68 @@
         self.deviatoricerror_resolution= 500.
         self.deviatoricerror_threshold = 0
         self.deviatoricerror_maximum    = 0
+        self.restart = 0.
         return self
     #}}}
     def checkconsistency(self,md,solution,analyses):    # {{{
@@ -83 +86 @@
         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 +107 @@
         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')
     # }}}

../trunk-jpl/src/m/classes/matestar.py

@@ +1 @@
+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')
+        # }}}

../trunk-jpl/src/m/classes/frictionsommers.py

@@ +1 @@
+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)     
+    # }}}

../trunk-jpl/src/m/classes/SMBgemb.py

@@ +1 @@
+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')
+        # }}}
+

../trunk-jpl/src/m/classes/taoinversion.py

@@ -5 +5 @@
 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 +99 @@
                 
                 #several responses can be used:
                 self.cost_functions=101;
-
                 return self
 
         def extrude(self,md):
@@ -118 +119 @@
                 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 +143 @@
                 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 +162 @@
                         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')

../trunk-jpl/src/m/classes/SMBgemb.m

@@ -357 +357 @@
                         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);

../trunk-jpl/src/m/classes/plumebasalforcings.py

@@ +1 @@
+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')
+        # }}}

../trunk-jpl/test/NightlyRun/test701.py

@@ +1 @@
+#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
+                ]

../trunk-jpl/test/NightlyRun/test703.py

@@ +1 @@
+#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,
+        ]

../trunk-jpl/test/NightlyRun/test293.py

@@ +1 @@
+#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,
+        ]

../trunk-jpl/src/m/mesh/bamgflowband.py

@@ +1 @@
+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

../trunk-jpl/src/m/mesh/bamg.py

@@ -1 +1 @@
 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 +79 @@
 
                 #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 +91 @@
                 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 +118 @@
                         #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 +331 @@
         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)