Index: issm/trunk-jpl/test/NightlyRun/test2004.m
===================================================================
--- issm/trunk-jpl/test/NightlyRun/test2004.m	(revision 24717)
+++ issm/trunk-jpl/test/NightlyRun/test2004.m	(revision 24717)
@@ -0,0 +1,1154 @@
+%Test Name: Earth_Antarctica_GIA
+
+%Data paths: {{{
+modeldatapath='/Users/larour/ModelData';
+shppath='/Users/larour/issm-jpl/proj-group/qgis/Slr';
+gshhsshapefile=[modeldatapath '/Gshhg/Shp/GSHHS_shp/c/GSHHS_c_L1.shp'];
+%}}}
+
+%create sealevel model to hold our information: 
+sl=sealevelmodel();
+
+%Create basins using boundaries from shapefile: %{{{
+sl.addbasin(basin('continent','southamerica','name','southamerica','epsg',5530,'boundaries',{... %{{{
+	boundary('shppath',shppath,'shpfilename','PanamaWest','epsg',4326),...
+	boundary('shppath',shppath,'shpfilename','SouthAmericaNorthWest','epsg',4326,'orientation','reverse'),...
+	boundary('shppath',shppath,'shpfilename','NorthAmericaSouthWest','epsg',4326),...
+	boundary('shppath',shppath,'shpfilename','SouthAmericaWest','epsg',4326),...
+	boundary('shppath',shppath,'shpfilename','MarieByrd1Summit','epsg',3031),...
+	boundary('shppath',shppath,'shpfilename','Ellsworth2','epsg',3031,'orientation','reverse'),...
+	boundary('shppath',shppath,'shpfilename','Ellsworth1Summit','epsg',3031),...
+	boundary('shppath',shppath,'shpfilename','Bellingshausen2','epsg',3031,'orientation','reverse'),...
+	boundary('shppath',shppath,'shpfilename','Bellingshausen1Summit','epsg',3031),...
+	boundary('shppath',shppath,'shpfilename','Wilkins1','epsg',3031),...
+	boundary('shppath',shppath,'shpfilename','Wilkins2Summit','epsg',3031),...
+	boundary('shppath',shppath,'shpfilename','Peninsula','epsg',3031),...
+	boundary('shppath',shppath,'shpfilename','Palmer1Summit','epsg',3031),...
+	boundary('shppath',shppath,'shpfilename','Palmer1','epsg',3031),...
+	boundary('shppath',shppath,'shpfilename','Ronne2Summit','epsg',3031),...
+	boundary('shppath',shppath,'shpfilename','Ronne2b','epsg',3031),...
+	boundary('shppath',shppath,'shpfilename','SlessorSummit','epsg',3031),...
+	boundary('shppath',shppath,'shpfilename','Slessor1','epsg',3031,'orientation','reverse'),...
+	boundary('shppath',shppath,'shpfilename','BruntSummit','epsg',3031),...
+	boundary('shppath',shppath,'shpfilename','Brunt1','epsg',3031),...
+	boundary('shppath',shppath,'shpfilename','DronningMaud1Summit','epsg',3031),...
+	boundary('shppath',shppath,'shpfilename','SouthAmericaEast','epsg',4326),...
+	boundary('shppath',shppath,'shpfilename','MidAtlantic','epsg',4326),...
+	boundary('shppath',shppath,'shpfilename','SouthAmericaNorthEast','epsg',4326,'orientation','reverse'),...
+	boundary('shppath',shppath,'shpfilename','PanamaEast','epsg',4326),...
+	boundary('shppath',shppath,'shpfilename','Panama','epsg',4326)}));
+%}}}
+sl.addbasin(basin('continent','northamerica','name','northamerica','epsg',3628,'boundaries',{... %{{{
+	boundary('shppath',shppath,'shpfilename','NorthAmericaSouthWest','epsg',4326),...
+	boundary('shppath',shppath,'shpfilename','SouthAmericaNorthWest','epsg',4326),...
+	boundary('shppath',shppath,'shpfilename','PanamaWest','epsg',4326),...
+	boundary('shppath',shppath,'shpfilename','Panama','epsg',4326,'orientation','reverse'),...
+	boundary('shppath',shppath,'shpfilename','PanamaEast','epsg',4326),...
+	boundary('shppath',shppath,'shpfilename','SouthAmericaNorthEast','epsg',4326),...
+	boundary('shppath',shppath,'shpfilename','MidAtlantic','epsg',4326),...
+	boundary('shppath',shppath,'shpfilename','Atlantic','epsg',4326),...
+	boundary('shppath',shppath,'shpfilename','South2Summit','epsg',3413),...
+	boundary('shppath',shppath,'shpfilename','South2','epsg',3413,'orientation','reverse'),...
+	boundary('shppath',shppath,'shpfilename','South1Summit','epsg',3413),...
+	boundary('shppath',shppath,'shpfilename','SouthWest1','epsg',3413,'orientation','reverse'),...
+	boundary('shppath',shppath,'shpfilename','SouthWest1Summit','epsg',3413),...
+	boundary('shppath',shppath,'shpfilename','Russell1','epsg',3413,'orientation','reverse'),...
+	boundary('shppath',shppath,'shpfilename','Russell1Summit','epsg',3413),...
+	boundary('shppath',shppath,'shpfilename','Jakobshavn1','epsg',3413,'orientation','reverse'),...
+	boundary('shppath',shppath,'shpfilename','Jakobshavn1Summit','epsg',3413),...
+	boundary('shppath',shppath,'shpfilename','Store1','epsg',3413,'orientation','reverse'),...
+	boundary('shppath',shppath,'shpfilename','Store1Summit','epsg',3413),...
+	boundary('shppath',shppath,'shpfilename','NorthWest1','epsg',3413,'orientation','reverse'),...
+	boundary('shppath',shppath,'shpfilename','NorthWest1Summit','epsg',3413),...
+	boundary('shppath',shppath,'shpfilename','Heilprin1','epsg',3413,'orientation','reverse'),...
+	boundary('shppath',shppath,'shpfilename','Petermann2Summit','epsg',3413),...
+	boundary('shppath',shppath,'shpfilename','Petermann1','epsg',3413,'orientation','reverse'),...
+	boundary('shppath',shppath,'shpfilename','Petermann1Summit','epsg',3413),...
+	boundary('shppath',shppath,'shpfilename','Ostenfeld1','epsg',3413,'orientation','reverse'),...
+	boundary('shppath',shppath,'shpfilename','Ostenfeld1Summit','epsg',3413),...
+	boundary('shppath',shppath,'shpfilename','Academy1','epsg',3413,'orientation','reverse'),...
+	boundary('shppath',shppath,'shpfilename','Academy1Summit','epsg',3413),...
+	boundary('shppath',shppath,'shpfilename','ArcticCanada','epsg',4326),...
+	boundary('shppath',shppath,'shpfilename','NorthAmericaNorthWest','epsg',4326),...
+	boundary('shppath',shppath,'shpfilename','NorthAmericaWest','epsg',4326)}));
+%}}}
+sl.addbasin(basin('continent','australia','name','australia','epsg',4462,'boundaries',{... %{{{
+	boundary('shppath',shppath,'shpfilename','Ross3Summit','epsg',3031),...
+	boundary('shppath',shppath,'shpfilename','Australia','epsg',4326,'orientation','reverse'),...
+	boundary('shppath',shppath,'shpfilename','SouthAsiaSummit','epsg',4326),...
+	boundary('shppath',shppath,'shpfilename','SouthAsia','epsg',4326,'orientation','reverse'),...
+	boundary('shppath',shppath,'shpfilename','IndianOceanSummit','epsg',4326),...
+	boundary('shppath',shppath,'shpfilename','IndianOcean','epsg',4326,'orientation','reverse'),...
+	boundary('shppath',shppath,'shpfilename','AmerySummit','epsg',3031),...
+	boundary('shppath',shppath,'shpfilename','QueenMary2','epsg',3031),...
+	boundary('shppath',shppath,'shpfilename','QueenMary1Summit','epsg',3031),...
+	boundary('shppath',shppath,'shpfilename','Wilkes3','epsg',3031,'orientation','reverse'),...
+	boundary('shppath',shppath,'shpfilename','Wilkes2Summit','epsg',3031),...
+	boundary('shppath',shppath,'shpfilename','Adelie2','epsg',3031,'orientation','reverse'),...
+	boundary('shppath',shppath,'shpfilename','Adelie1Summit','epsg',3031),...
+	boundary('shppath',shppath,'shpfilename','Victoria2','epsg',3031,'orientation','reverse'),...
+	boundary('shppath',shppath,'shpfilename','Victoria1Summit','epsg',3031),...
+	boundary('shppath',shppath,'shpfilename','Ross4','epsg',3031,'orientation','reverse')}));
+%}}}
+sl.addbasin(basin('continent','eurasia','name','eurasia','epsg',3416,'boundaries',{... %{{{
+	boundary('shppath',shppath,'shpfilename','South2Summit','epsg',3413),...
+	boundary('shppath',shppath,'shpfilename','Atlantic','epsg',4326,'orientation','reverse'),...
+	boundary('shppath',shppath,'shpfilename','MidAtlantic','epsg',4326),...
+	boundary('shppath',shppath,'shpfilename','SouthAmericaEast','epsg',4326,'orientation','reverse'),...
+	boundary('shppath',shppath,'shpfilename','DronningMaud1Summit','epsg',3031),...
+	boundary('shppath',shppath,'shpfilename','DronningMaud1','epsg',3031,'orientation','reverse'),...
+	boundary('shppath',shppath,'shpfilename','QueenMaud2Summit','epsg',3031),...
+	boundary('shppath',shppath,'shpfilename','QueenMaud2','epsg',3031,'orientation','reverse'),...
+	boundary('shppath',shppath,'shpfilename','Kemp2Summit','epsg',3031),...
+	boundary('shppath',shppath,'shpfilename','Kemp2','epsg',3031,'orientation','reverse'),...
+	boundary('shppath',shppath,'shpfilename','KempSummit','epsg',3031),...
+	boundary('shppath',shppath,'shpfilename','Amery5','epsg',3031),...
+	boundary('shppath',shppath,'shpfilename','AmerySummit','epsg',3031),...
+	boundary('shppath',shppath,'shpfilename','IndianOcean','epsg',4326),...
+	boundary('shppath',shppath,'shpfilename','IndianOceanSummit','epsg',4326),...
+	boundary('shppath',shppath,'shpfilename','SouthAsia','epsg',4326),...
+	boundary('shppath',shppath,'shpfilename','SouthAsiaSummit','epsg',4326),...
+	boundary('shppath',shppath,'shpfilename','ArcticRussia','epsg',4326,'orientation','reverse'),...
+	boundary('shppath',shppath,'shpfilename','NorthAmericaNorthWest','epsg',4326),...
+	boundary('shppath',shppath,'shpfilename','ArcticCanada','epsg',4326,'orientation','reverse'),...
+	boundary('shppath',shppath,'shpfilename','Academy1Summit','epsg',3413),...
+	boundary('shppath',shppath,'shpfilename','NEGIS1','epsg',3413,'orientation','reverse'),...
+	boundary('shppath',shppath,'shpfilename','NEGIS1Summit','epsg',3413),...
+	boundary('shppath',shppath,'shpfilename','Bistrup2','epsg',3413,'orientation','reverse'),...
+	boundary('shppath',shppath,'shpfilename','Bistrup1Summit','epsg',3413),...
+	boundary('shppath',shppath,'shpfilename','Daugaard-Jensen5','epsg',3413,'orientation','reverse'),...
+	boundary('shppath',shppath,'shpfilename','Daugaard-Jensen-Summit4','epsg',3413),...
+	boundary('shppath',shppath,'shpfilename','Geikie2','epsg',3413,'orientation','reverse'),...
+	boundary('shppath',shppath,'shpfilename','Geikie1Summit','epsg',3413),...
+	boundary('shppath',shppath,'shpfilename','Kangerlussuaq3','epsg',3413,'orientation','reverse'),...
+	boundary('shppath',shppath,'shpfilename','Kangerlussuaq2Summit','epsg',3413),...
+	boundary('shppath',shppath,'shpfilename','Helheim3','epsg',3413,'orientation','reverse'),...
+	boundary('shppath',shppath,'shpfilename','Helheim2Summit','epsg',3413),...
+	boundary('shppath',shppath,'shpfilename','KogeBugt3','epsg',3413,'orientation','reverse'),...
+	boundary('shppath',shppath,'shpfilename','KogeBugt2Summit','epsg',3413),...
+	boundary('shppath',shppath,'shpfilename','SouthEast3','epsg',3413,'orientation','reverse'),...
+	boundary('shppath',shppath,'shpfilename','SouthEast2Summit','epsg',3413),...
+	boundary('shppath',shppath,'shpfilename','South3','epsg',3413,'orientation','reverse')}));
+%}}}
+sl.addbasin(basin('continent','pacific','name','pacific','epsg',3031,'boundaries',{... %{{{
+	boundary('shppath',shppath,'shpfilename','Ross3Summit','epsg',3031),...
+	boundary('shppath',shppath,'shpfilename','MarieByrd2','epsg',3031,'orientation','reverse'),...
+	boundary('shppath',shppath,'shpfilename','MarieByrd1Summit','epsg',3031),...
+	boundary('shppath',shppath,'shpfilename','SouthAmericaWest','epsg',4326,'orientation','reverse'),...
+	boundary('shppath',shppath,'shpfilename','NorthAmericaSouthWest','epsg',4326),...
+	boundary('shppath',shppath,'shpfilename','NorthAmericaWest','epsg',4326,'orientation','reverse'),...
+	boundary('shppath',shppath,'shpfilename','NorthAmericaNorthWest','epsg',4326),...
+	boundary('shppath',shppath,'shpfilename','ArcticRussia','epsg',4326),...
+	boundary('shppath',shppath,'shpfilename','SouthAsiaSummit','epsg',4326),...
+	boundary('shppath',shppath,'shpfilename','Australia','epsg',4326)}));
+%}}}
+sl.addbasin(basin('continent','greenland','name','greenland','epsg',3413,'boundaries',{... %{{{
+	boundary('shppath',shppath,'shpfilename','South2Summit','epsg',3413),...
+	boundary('shppath',shppath,'shpfilename','South2','epsg',3413,'orientation','reverse'),...
+	boundary('shppath',shppath,'shpfilename','South1Summit','epsg',3413),...
+	boundary('shppath',shppath,'shpfilename','SouthWest1','epsg',3413,'orientation','reverse'),...
+	boundary('shppath',shppath,'shpfilename','SouthWest1Summit','epsg',3413),...
+	boundary('shppath',shppath,'shpfilename','Russell1','epsg',3413,'orientation','reverse'),...
+	boundary('shppath',shppath,'shpfilename','Russell1Summit','epsg',3413),...
+	boundary('shppath',shppath,'shpfilename','Jakobshavn1','epsg',3413,'orientation','reverse'),...
+	boundary('shppath',shppath,'shpfilename','Jakobshavn1Summit','epsg',3413),...
+	boundary('shppath',shppath,'shpfilename','Store1','epsg',3413,'orientation','reverse'),...
+	boundary('shppath',shppath,'shpfilename','Store1Summit','epsg',3413),...
+	boundary('shppath',shppath,'shpfilename','NorthWest1','epsg',3413,'orientation','reverse'),...,
+	boundary('shppath',shppath,'shpfilename','NorthWest1Summit','epsg',3413),...
+	boundary('shppath',shppath,'shpfilename','Heilprin1','epsg',3413,'orientation','reverse'),...,
+	boundary('shppath',shppath,'shpfilename','Petermann2Summit','epsg',3413),...
+	boundary('shppath',shppath,'shpfilename','Petermann1','epsg',3413,'orientation','reverse'),...,
+	boundary('shppath',shppath,'shpfilename','Petermann1Summit','epsg',3413),...
+	boundary('shppath',shppath,'shpfilename','Ostenfeld1','epsg',3413,'orientation','reverse'),...,
+	boundary('shppath',shppath,'shpfilename','Ostenfeld1Summit','epsg',3413),...
+	boundary('shppath',shppath,'shpfilename','Academy1','epsg',3413,'orientation','reverse'),...,
+	boundary('shppath',shppath,'shpfilename','Academy1Summit','epsg',3413),...
+	boundary('shppath',shppath,'shpfilename','NEGIS1','epsg',3413,'orientation','reverse'),...,
+	boundary('shppath',shppath,'shpfilename','NEGIS1Summit','epsg',3413),...
+	boundary('shppath',shppath,'shpfilename','Bistrup2','epsg',3413,'orientation','reverse'),...,
+	boundary('shppath',shppath,'shpfilename','Bistrup1Summit','epsg',3413),...
+	boundary('shppath',shppath,'shpfilename','Daugaard-Jensen5','epsg',3413,'orientation','reverse'),...,
+	boundary('shppath',shppath,'shpfilename','Daugaard-Jensen-Summit4','epsg',3413),...
+	boundary('shppath',shppath,'shpfilename','Geikie2','epsg',3413,'orientation','reverse'),...,
+	boundary('shppath',shppath,'shpfilename','Geikie1Summit','epsg',3413),...
+	boundary('shppath',shppath,'shpfilename','Kangerlussuaq3','epsg',3413,'orientation','reverse'),...,
+	boundary('shppath',shppath,'shpfilename','Kangerlussuaq2Summit','epsg',3413),...
+	boundary('shppath',shppath,'shpfilename','Helheim3','epsg',3413,'orientation','reverse'),...,
+	boundary('shppath',shppath,'shpfilename','Helheim2Summit','epsg',3413),...
+	boundary('shppath',shppath,'shpfilename','KogeBugt3','epsg',3413,'orientation','reverse'),...,
+	boundary('shppath',shppath,'shpfilename','KogeBugt2Summit','epsg',3413),...
+	boundary('shppath',shppath,'shpfilename','SouthEast3','epsg',3413,'orientation','reverse'),...,
+	boundary('shppath',shppath,'shpfilename','SouthEast2Summit','epsg',3413),...
+	boundary('shppath',shppath,'shpfilename','South3','epsg',3413,'orientation','reverse')}));
+%}}}
+sl.addbasin(basin('continent','antarctica','name','antarctica','epsg',3031,'boundaries',{... %{{{
+		boundary('shppath',shppath,'shpfilename','DronningMaud1Summit','epsg',3031),...
+		boundary('shppath',shppath,'shpfilename','DronningMaud1','epsg',3031,'orientation','reverse'),...
+		boundary('shppath',shppath,'shpfilename','QueenMaud2Summit','epsg',3031),...
+		boundary('shppath',shppath,'shpfilename','QueenMaud2','epsg',3031,'orientation','reverse'),...
+		boundary('shppath',shppath,'shpfilename','Kemp2Summit','epsg',3031),...
+		boundary('shppath',shppath,'shpfilename','Kemp2','epsg',3031,'orientation','reverse'),...
+		boundary('shppath',shppath,'shpfilename','KempSummit','epsg',3031),...
+		boundary('shppath',shppath,'shpfilename','Amery5','epsg',3031),...
+		boundary('shppath',shppath,'shpfilename','AmerySummit','epsg',3031),...
+		boundary('shppath',shppath,'shpfilename','QueenMary2','epsg',3031),...
+		boundary('shppath',shppath,'shpfilename','QueenMary1Summit','epsg',3031),...
+		boundary('shppath',shppath,'shpfilename','Wilkes3','epsg',3031,'orientation','reverse'),...
+		boundary('shppath',shppath,'shpfilename','Wilkes2Summit','epsg',3031),...
+		boundary('shppath',shppath,'shpfilename','Adelie2','epsg',3031,'orientation','reverse'),...
+		boundary('shppath',shppath,'shpfilename','Adelie1Summit','epsg',3031),...
+		boundary('shppath',shppath,'shpfilename','Victoria2','epsg',3031,'orientation','reverse'),...
+		boundary('shppath',shppath,'shpfilename','Victoria1Summit','epsg',3031),...
+		boundary('shppath',shppath,'shpfilename','Ross4','epsg',3031,'orientation','reverse'),...
+		boundary('shppath',shppath,'shpfilename','Ross3Summit','epsg',3031),...
+		boundary('shppath',shppath,'shpfilename','MarieByrd2','epsg',3031,'orientation','reverse'),...
+		boundary('shppath',shppath,'shpfilename','MarieByrd1Summit','epsg',3031),...
+		boundary('shppath',shppath,'shpfilename','Ellsworth2','epsg',3031,'orientation','reverse'),...
+		boundary('shppath',shppath,'shpfilename','Ellsworth1Summit','epsg',3031),...
+		boundary('shppath',shppath,'shpfilename','Bellingshausen2','epsg',3031,'orientation','reverse'),...
+		boundary('shppath',shppath,'shpfilename','Bellingshausen1Summit','epsg',3031),...
+		boundary('shppath',shppath,'shpfilename','Wilkins1','epsg',3031),...
+		boundary('shppath',shppath,'shpfilename','Wilkins2Summit','epsg',3031),...
+		boundary('shppath',shppath,'shpfilename','Peninsula','epsg',3031),...
+		boundary('shppath',shppath,'shpfilename','Palmer1Summit','epsg',3031),...
+		boundary('shppath',shppath,'shpfilename','Palmer1','epsg',3031),...
+		boundary('shppath',shppath,'shpfilename','Ronne2Summit','epsg',3031),...
+		boundary('shppath',shppath,'shpfilename','Ronne2b','epsg',3031),...
+		boundary('shppath',shppath,'shpfilename','SlessorSummit','epsg',3031),...
+		boundary('shppath',shppath,'shpfilename','Slessor1','epsg',3031,'orientation','reverse'),...
+		boundary('shppath',shppath,'shpfilename','BruntSummit','epsg',3031),...
+		boundary('shppath',shppath,'shpfilename','Brunt1','epsg',3031)}));
+	%}}}
+%}}}
+%Go through basins and mesh:  %{{{
+%meshing parameters:  {{{
+hmin=500; hmax=2000; hmin=hmin*1000; hmax=hmax*1000; 
+tolerance=100; %tolerance of 100m on Earth position when mergin 3d meshes
+threshold=1;
+defaultoptions={'KeepVertices',0,'MaxCornerAngle',0.0000000001,'NoBoundaryRefinment',1}; 
+alreadyloaded=0;
+%}}}
+for i=sl.basinindx('basin','all'),
+
+	bas=sl.basins{i};
+	disp(sprintf('Meshing basin %s\n',bas.name));
+
+	%recover basin domain: 
+	domain=bas.contour();
+
+	%recover coastline inside basin, using GSHHS_c_L1. It's a lat,long file, hence epsg 4326
+	coastline=bas.shapefilecrop('shapefile',gshhsshapefile,'epsgshapefile',4326,'threshold',threshold); 
+
+	%mesh: 
+	md=bamg(model,'domain',domain,'subdomains',coastline,'hmin',hmin,'hmax',hmax,defaultoptions{:});
+	plotmodel(md,'data','mesh');pause(1);
+	error;
+
+	%miscellaneous: 
+	md.mesh.epsg=bas.epsg; md.miscellaneous.name=bas.name;
+
+	%recover mask where we have land: 
+	md.private.bamg.landmask=double(md.private.bamg.mesh.Triangles(:,4)>=1);
+
+	%vertex connectivity:
+	md.mesh.vertexconnectivity=NodeConnectivity(md.mesh.elements,md.mesh.numberofvertices);
+
+	%add model to sl icecaps:
+	sl.addicecap(md);
+end
+%}}}
+error;
+%Parameterize ice sheets : {{{
+
+for ind=sl.basinindx('continent',{'greenland','antarctica'}),
+	disp(sprintf('Parameterizing basin %s\n', sl.icecaps{ind}.miscellaneous.name));
+
+	md=sl.icecaps{ind}; bas=sl.basins{ind}; 
+
+	%basin/continent specific code: {{{
+	if bas.iscontinentany('antarctica'), 
+		openstreet=openstreetantarctica; 
+		domainoutline=domainoutlinea;
+	end; 
+	if bas.iscontinentany('greenland'), 
+		openstreet=openstreetgreenland; 
+		domainoutline=domainoutlineg;
+	end
+	%}}}
+	%masks :  %{{{
+	%new type of mask: 
+	md.mask=maskpsl; 
+
+	%land and ocean from open street: 
+	flags=1-(~ContourToNodes(md.mesh.x,md.mesh.y,openstreet,0));
+	pos=find(flags==0); flags(pos)=-1;
+	md.mask.land_levelset=flags;
+	md.mask.ocean_levelset=-flags;
+
+	%ice levelset from domain outlines: 
+	flags= -(~ContourToNodes(md.mesh.x,md.mesh.y,[shppath '/' domainoutline],1));
+	mds=extract(md,[shppath '/' domainoutline]);
+	flags(mds.mesh.extractedvertices)=~mds.mesh.vertexonboundary;
+	pos=find(flags==0 & md.mesh.vertexonboundary); flags(pos)=1;
+	md.mask.ice_levelset=-flags;
+
+	%on ice front, we are not on ocean; 
+	pos=find(md.mask.ice_levelset==0); md.mask.ocean_levelset(pos)=-1; md.mask.land_levelset(pos)=1;
+
+	%now, glaciers: 
+	mesh_glacier_levelset=md.private.bamg.landmask; %see meshing for why
+	md.mask.glacier_levelset=zeros(md.mesh.numberofvertices,1);
+	pos=find(mesh_glacier_levelset); 
+	md.mask.glacier_levelset(md.mesh.elements(pos,:))=1;
+
+	pos=find(md.mask.glacier_levelset); 
+	if ~isempty(pos),
+		mds=extract(md,md.mask.glacier_levelset);
+		md.mask.ice_levelset(mds.mesh.extractedvertices)=-(~mds.mesh.vertexonboundary);
+	end
+	%}}} 
+	%initial grounding line position: {{{
+	if bas.iscontinentany('antarctica'), % {{{
+
+		%figure out initial grounding line position from the bedmap2 dataset: 
+		land_type=interpBedmap2(md.mesh.x,md.mesh.y,'icemask_grounded_and_shelves');
+
+		pos=find(isnan(land_type)); 
+		pos2=find(~isnan(land_type)); 
+		for i=1:length(pos),
+			in=find_point(md.mesh.x(pos2),md.mesh.y(pos2),md.mesh.x(pos(i)),md.mesh.y(pos(i)));
+			land_type(pos(i))=land_type(pos2(in));
+		end
+
+		gridoniceshelf=zeros(md.mesh.numberofvertices,1);
+		gridoniceshelf(land_type>0.9)=land_type(land_type>0.9);
+		gridoniceshelf(gridoniceshelf>0)=-1;
+		gridoniceshelf(gridoniceshelf~=-1)=1;
+
+		md.mask.groundedice_levelset=gridoniceshelf;
+
+		%correction to land and ocean levelset: ice shelf is not on land! 
+		pos=find(md.mask.ice_levelset<=0 & md.mask.groundedice_levelset<=0);
+		md.mask.ocean_levelset(pos)=1;
+		md.mask.land_levelset(pos)=-1;
+		%}}}
+	elseif bas.iscontinentany('greenland'), % {{{
+
+		mask=interpBedmachine(md.mesh.x,md.mesh.y,'mask','greenland');
+
+		gridoniceshelf=ones(md.mesh.numberofvertices,1);
+		pos=find(md.mask.ice_levelset>0); groundedice_levelset(pos)=-1;
+		pos=find(mask<=0); gridoniceshelf(pos)=0;
+
+		md.mask.groundedice_levelset=gridoniceshelf; 
+		%}}}
+	end
+	%}}}
+	%latlong:  % {{{
+	[md.mesh.long,md.mesh.lat]=gdaltransform(md.mesh.x,md.mesh.y,sprintf('EPSG:%i',md.mesh.epsg),'EPSG:4326'); 
+	%}}}
+	%flowequation: % {{{
+	md=setflowequation(md,'SSA','all');
+	%}}}
+	%geometry: {{{
+	if bas.iscontinentany('antarctica'),% {{{
+		di=md.materials.rho_ice/md.materials.rho_water;
+
+		disp('      reading dem');
+		surface=interpBedmachine(md.mesh.x,md.mesh.y,'surface','antarctica');
+		surface=griddata(md.mesh.x(~isnan(surface)),md.mesh.y(~isnan(surface)),surface(~isnan(surface)),md.mesh.x,md.mesh.y,'nearest');
+		surfaceold=surface;
+
+		disp('      reading firn layer');
+		load(firnpath)
+		firn_layer=InterpFromGridToMesh(x1,y1,firn,md.mesh.x,md.mesh.y,0);
+
+		disp('      reading bedrock');
+		base=interpBedmachine(md.mesh.x,md.mesh.y,'bed','antarctica');
+		bedmap=interpBedmap2(md.mesh.x,md.mesh.y,'bed');
+		base(base<-8000 | isnan(base))=bedmap(base<-8000 | isnan(base));
+
+		%thickness: 
+		thickness=surface-base;
+
+		%hydrostatic equilibrium on ice shelves:
+		pos=find(md.mask.groundedice_levelset<0); thickness(pos)=surface(pos)/(1-di)-firn_layer(pos);
+		surface(pos)=(1-di)*thickness(pos);
+		base(pos)=-di*thickness(pos);
+
+		%firn layer over grounded ice: 
+		pos=find(md.mask.groundedice_levelset>=0); thickness(pos)=thickness(pos)-firn_layer(pos);
+
+		%check again: 
+		pos=find(thickness<1); thickness(pos)=1;
+
+		%reset: 
+		surface=base+thickness;
+
+		%make sure we are at hydrostatic at the grounding line: 
+		bed=base;bed(md.mask.groundedice_levelset<0)=base(md.mask.groundedice_levelset<0)-2000;
+
+		%set: 
+		md.geometry.surface=surface;
+		md.geometry.bed=bed;
+		md.geometry.base=base;
+		md.geometry.thickness=thickness;
+		%}}}
+	elseif bas.iscontinentany('greenland'),% {{{
+		mask=interpBedmachine(md.mesh.x,md.mesh.y,'mask','greenland');
+		surface=interpBedmachine(md.mesh.x,md.mesh.y,'surface','greenland');
+		thickness=interpBedmachine(md.mesh.x,md.mesh.y,'thickness','greenland');
+		base=surface-thickness;
+		di=md.materials.rho_ice/md.materials.rho_water;
+
+		pos=find(isnan(base) | isnan(surface));
+		thickness(pos)=1; 
+		surface(pos)=(1-di)*thickness(pos); 
+		base(pos)=-di*thickness(pos);
+
+		%places where thickness is 0: 
+		pos=find(thickness==0);
+		thickness(pos)=.1; surface(pos)=base(pos)+thickness(pos);
+
+		%adjust bed and base: 
+		bed=base;
+		pos=find(md.mask.groundedice_levelset>=0); bed(pos)=base(pos);
+		pos=find(md.mask.groundedice_levelset<0); bed(pos)=base(pos)-100;
+
+		md.geometry.surface=surface;
+		md.geometry.bed=bed;
+		md.geometry.base=bed;
+		md.geometry.thickness=thickness;
+
+	end %}}}
+	% }}}
+	%velocities:  %{{{
+	if bas.iscontinentany('antarctica'),% {{{
+
+		velnc        =[jplsvn '/database/antarctica/velocity_ref_v28Apr2011bamber_900m_median.nc'];
+
+		%read and process netcdf
+		xmin    = double(ncreadatt(velnc,'/','xmin'));
+		ymax    = double(ncreadatt(velnc,'/','ymax'));
+		nx      = double(ncreadatt(velnc,'/','nx'));
+		ny      = double(ncreadatt(velnc,'/','ny'));
+		spacing = double(ncreadatt(velnc,'/','spacing'));
+		vx      = double(ncread(velnc,'vx'));
+		vy      = double(ncread(velnc,'vy'));
+		x=xmin+(0:1:nx)'*spacing;
+		y=(ymax-ny*spacing)+(0:1:ny)'*spacing;
+		vx=InterpFromGridToMesh(x,y,flipud(vx'),md.mesh.x,md.mesh.y,0);
+		vy=InterpFromGridToMesh(x,y,flipud(vy'),md.mesh.x,md.mesh.y,0);
+		vel=sqrt(vx.^2+vy.^2);
+
+		%}}}
+	elseif bas.iscontinentany('greenland'),% {{{
+
+		velocities=refinementmetric(md.mesh,'greenland','greenland','nsidcvxvy'); 
+		vx=velocities(:,1); vy=velocities(:,2); vel=sqrt(vx.^2+vy.^2);
+
+		velocities=refinementmetric(md.mesh,'greenland','greenland','joughinvxvy'); 
+		vxfar=velocities(:,1); vyfar=velocities(:,2); velfar=sqrt(vxfar.^2+vyfar.^2);
+
+		%use velfar inland (from mosaic), and 2008 velocities nearcoastline: 
+		pos=find(vel==0 | isnan(vel)); vel(pos)=velfar(pos); vx(pos)=vxfar(pos);vy(pos)=vyfar(pos);
+
+		%set water velocity to 0: 
+		pos=find(isnan(vel)); vel(pos)=0; vx(pos)=0; vy(pos)=0;
+
+	end % }}}
+
+	md.inversion.vx_obs=vx;
+	md.inversion.vy_obs=vy;
+	md.inversion.vel_obs=vel;
+	md.initialization.vx=md.inversion.vx_obs;
+	md.initialization.vy=md.inversion.vy_obs;
+	md.initialization.vz=zeros(md.mesh.numberofvertices,1);
+	md.initialization.vel=md.inversion.vel_obs;
+
+	% }}} 
+	%friction:  {{{
+	[sx,sy,s]=slope(md); sslope=averaging(md,s,0);
+
+	pos=find(md.inversion.vel_obs==0);
+	vel=max(md.inversion.vel_obs,0.1); vel(pos)=1;
+
+	md.friction.coefficient=sqrt(md.materials.rho_ice*md.geometry.thickness.*(sslope)./((md.materials.rho_ice*md.geometry.thickness+md.materials.rho_water*md.geometry.base).*vel/md.constants.yts));
+	md.friction.coefficient=min(md.friction.coefficient,200);
+
+	min_drag_coeff=35;
+	min_drag_coeff_outlets=5;
+	threshhold_drag_coeff=70;
+	background_drag_coeff=200;
+
+	maxvel=30;
+	pos=find(md.friction.coefficient<min_drag_coeff);
+	md.friction.coefficient(pos)=min_drag_coeff;
+
+	%Take care of iceshelves: no drag md.drag
+	pos=find(md.mask.groundedice_levelset<0);
+	md.friction.coefficient(pos)=0;
+	md.friction.p=ones(md.mesh.numberofelements,1);
+	md.friction.q=ones(md.mesh.numberofelements,1);
+	% }}}
+	%Temperatures and surface mass balance: {{{ 
+	if bas.iscontinentany('antarctica'),% {{{
+		%smb {{{
+		smbpath = [modeldatapath2 '/RACMO2Antarctica/'];
+
+		Years = 1979:2010; %we are going to need 1980 -> 1990 as average, and perturbation from 2005 -> 2010
+
+		ncdata = sprintf('%s/%s-%i%s',smbpath,'smb.KNMI',1979,'.ANT3K27.DRIFTALB.nc');
+		lat = ncread(ncdata,'lat'); lon = ncread(ncdata,'lon');
+		lat = double(lat); lon = double(lon);
+		lat = lat(:); lon = lon(:); 
+		[x,y]=ll2xy(lat,lon,-1);
+		index=BamgTriangulate(x,y);
+		di=md.materials.rho_ice/md.materials.rho_water;
+
+		smbs=zeros(length(x),numel(Years));
+		for i=1:length(Years),
+			year=Years(i);
+			ncdata = sprintf('%s/%s-%i%s',smbpath,'smb.KNMI',year,'.ANT3K27.DRIFTALB.nc');
+			smb=ncread(ncdata,'smb'); smb=sum(smb(:,:,:),3)/1000/di;  smb=smb(:);
+			smbs(:,i)=smb;
+		end
+		smb=InterpFromMeshToMesh2d(index,x,y,smbs,md.mesh.x,md.mesh.y);
+
+		%mean from 1979 to 1990: 
+		pos=find(Years>=1979 & Years<=1990);
+		smbmean7990=mean(smb(:,pos),2);
+		smbmean7990=griddata(md.mesh.x(~isnan(smbmean7990)),md.mesh.y(~isnan(smbmean7990)),smbmean7990(~isnan(smbmean7990)),md.mesh.x,md.mesh.y,'nearest');
+
+		%mean from  2005-2010:
+		pos=find(Years>=2005 & Years<=2010);
+		smbmean0510=mean(smb(:,pos),2);
+		smbmean0510=griddata(md.mesh.x(~isnan(smbmean0510)),md.mesh.y(~isnan(smbmean0510)),smbmean0510(~isnan(smbmean0510)),md.mesh.x,md.mesh.y,'nearest');
+
+		md.smb.mass_balance = smbmean0510;;
+
+		pos=find(md.mask.glacier_levelset);
+		md.smb.mass_balance(pos)=smbmean0510(pos)-smbmean7990(pos);
+
+		%}}}
+		%temperature:  {{{
+		temperaturepath=['/Users/larour/ModelData/RACMO2Antarctica/annt33.mat'];
+		load(temperaturepath);
+		index=delaunay(x1,y1);
+		md.initialization.temperature=InterpFromMeshToMesh2d(index,x1(:),y1(:),annt33(:),md.mesh.x,md.mesh.y,'default',273.15);
+		md.initialization.pressure=md.materials.rho_ice*md.constants.g*(md.geometry.surface-md.geometry.base);
+		% }}}
+		%}}}
+	elseif bas.iscontinentany('greenland'),% {{{
+		%% Annual mean for the period given by Years (smb0) 
+		Years = 1960:2014; %we are going to need 1960 -> 1990 as average, and perturbation from 2010 -> 2014.
+
+		smb2 = nan(md.mesh.numberofvertices,numel(Years));
+		st2  = nan(md.mesh.numberofvertices,numel(Years));
+
+		%triangulate first: 
+		ncdata = [MARname '/MARv3.5-15km-monthly-ERA-40-1971.nc'];
+		lat = ncread(ncdata,'LAT'); lon = ncread(ncdata,'LON');
+		lat = double(lat); lon = double(lon);
+		lat = lat(:); lon = lon(:); 
+		[x,y]=ll2xy(lat,lon,+1,45,70);
+		index=BamgTriangulate(x,y);
+
+		smb1s=zeros(length(x),numel(Years));
+		st1s=zeros(length(x),numel(Years));
+
+		for ii = 1:numel(Years)
+			if Years(ii)<=1978
+				ncdata = [MARname '/MARv3.5-15km-monthly-ERA-40-' num2str(Years(ii)) '.nc'];
+			elseif Years(ii)>1978
+				ncdata = [MARname '/MARv3.5-15km-monthly-ERA-Interim-' num2str(Years(ii)) '.nc'];
+			end
+
+			st1 = ncread(ncdata,'STcorr'); % Surface temperature 
+			st1(st1<-0.9000e30)=nan;
+			st1 = double(mean(st1,3)); % Annual mean surface temp
+
+			st1s(:,ii)=st1(:);
+
+			smb1 = ncread(ncdata,'SMBcorr'); % in mmWe/month
+			smb1(smb1<-0.900e30)=nan;
+			smb1 = double(sum(smb1,3)/1000); %mWe/year
+
+			smb1s(:,ii)=smb1(:);
+
+		end
+
+		st2s=InterpFromMeshToMesh2d(index,x,y,st1s,md.mesh.x,md.mesh.y);         
+		smb2s=InterpFromMeshToMesh2d(index,x,y,smb1s,md.mesh.x,md.mesh.y);
+
+		smb19601990=mean(smb2s(:,1:31),2);
+		smb20102014=mean(smb2s(:,end-4:end),2);
+
+		st20102014 = mean(st2s(:,end-4:end),2); 
+
+		md.initialization.temperature=273.25+st20102014; 
+		md.initialization.pressure=md.materials.rho_ice*md.constants.g*(md.geometry.surface-md.geometry.base);
+		md.smb.mass_balance = smb20102014;
+
+		pos=find(md.mask.glacier_levelset);
+		md.smb.mass_balance(pos)=smb20102014(pos)-smb19601990(pos);
+
+	end % }}}
+	%}}}
+	%Mechanical boundary conditions: {{{
+	%intialize spc boundary conditions:
+	md.stressbalance.spcvx=NaN*ones(md.mesh.numberofvertices,1);
+	md.stressbalance.spcvy=NaN*ones(md.mesh.numberofvertices,1);
+	md.stressbalance.spcvz=NaN*ones(md.mesh.numberofvertices,1);
+	md.stressbalance.referential=NaN*ones(md.mesh.numberofvertices,6);
+
+	%boundary conditions according to masks: 
+
+	%water (freeze velocity to 0)
+	pos=find(md.mask.ocean_levelset>=0);
+	md.stressbalance.spcvx(pos)=0;
+	md.stressbalance.spcvy(pos)=0;
+	md.stressbalance.spcvz(pos)=0;
+
+	%glaciers (freeeze velocity to observed)
+	pos=find(md.mask.glacier_levelset);
+	md.stressbalance.spcvx(pos)=0; %flux divergence will be equated to smb perturbation.
+	md.stressbalance.spcvy(pos)=0;
+
+	md.stressbalance.spcvz(pos)=0;
+
+	%land with no ice (freeeze velocity to 0)
+	pos=find(md.mask.land_levelset>=0 & md.mask.ice_levelset>0 & ~md.mask.glacier_levelset);
+	md.stressbalance.spcvx(pos)=0;
+	md.stressbalance.spcvy(pos)=0;
+	md.stressbalance.spcvz(pos)=0;
+
+	%if strictly ice (and not glacier), then no boundary condition: 
+	pos=find(md.mask.ice_levelset<=0 & ~md.mask.glacier_levelset); 
+	md.stressbalance.spcvx(pos)=NaN;
+	md.stressbalance.spcvy(pos)=NaN;
+	md.stressbalance.spcvz(pos)=NaN;
+
+	%constrain boundaries of basin to observed velocities: 
+	pos=find(md.mesh.vertexonboundary & md.mask.ice_levelset<=0); 
+	md.stressbalance.spcvx(pos)=md.inversion.vx_obs(pos);
+	md.stressbalance.spcvy(pos)=md.inversion.vy_obs(pos);
+	md.stressbalance.spcvz(pos)=0;
+
+	%fixing quirks here and there: 
+	if strcmpi(bas.name,'Ross'),
+		flags = ContourToNodes(mdband.mesh.x,mdband.mesh.y,[antarcticashppath '/Murdo.exp'],1);
+		pos=find(flags); 
+
+		md.stressbalance.spcvx(pos)=md.inversion.vx_obs(pos);
+		md.stressbalance.spcvy(pos)=md.inversion.vy_obs(pos);
+		md.stressbalance.spcvz(pos)=0;
+	end
+
+	%prognostic
+	md.masstransport.spcthickness=NaN*ones(md.mesh.numberofvertices,1);
+	% }}}
+	%Thermal boundary conditions: {{{
+	if bas.iscontinentany('antarctica'),% {{{
+		searisepath  =[modeldatapath '/SeaRISE/Antarctica5km_shelves_v1.0'];
+		heatfluxpath   =[searisepath '/bheatflx_fox.mat'];
+
+		if strcmpi(bas.name,'ellsworth'), % {{{
+			md.basalforcings=linearbasalforcings(md);
+			md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
+			md.basalforcings.deepwater_melting_rate=79.3397;
+			md.basalforcings.deepwater_elevation=-1000;
+			md.basalforcings.upperwater_elevation=-141.59; 
+			%}}}
+		elseif strcmpi(bas.name,'amery'), % {{{
+			md.basalforcings=linearbasalforcings(md);
+			md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
+			md.basalforcings.deepwater_melting_rate=8;
+			md.basalforcings.deepwater_elevation=-1500;
+			md.basalforcings.upperwater_elevation=-400; % }}}
+		elseif strcmpi(bas.name,'ross'), % {{{
+			md.basalforcings=linearbasalforcings(md);
+			md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
+			md.basalforcings.deepwater_melting_rate=5;
+			md.basalforcings.deepwater_elevation=-1100;
+			md.basalforcings.upperwater_elevation=-200; % }}}
+		elseif strcmpi(bas.name,'ronne'), % {{{
+			md.basalforcings=linearbasalforcings(md);
+			md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
+			md.basalforcings.deepwater_melting_rate=4;
+			md.basalforcings.deepwater_elevation=-1800;
+			md.basalforcings.upperwater_elevation=-400; % }}}
+		elseif strcmpi(bas.name,'brunt'), % {{{
+			md.basalforcings=linearbasalforcings(md);
+			md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
+			md.basalforcings.deepwater_melting_rate=.3;
+			md.basalforcings.deepwater_elevation=-900;
+			md.basalforcings.upperwater_elevation=-100; 
+			% }}}
+		elseif strcmpi(bas.name,'kemp'), % {{{
+			md.basalforcings=linearbasalforcings(md);
+			md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
+			md.basalforcings.deepwater_melting_rate=.2;
+			md.basalforcings.deepwater_elevation=-800;
+			md.basalforcings.upperwater_elevation=-100; % }}}
+		elseif strcmpi(bas.name,'queenMaud'), % {{{
+			md.basalforcings=linearbasalforcings(md);
+			md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
+			md.basalforcings.deepwater_melting_rate=.5;
+			md.basalforcings.deepwater_elevation=-650;
+			md.basalforcings.upperwater_elevation=-180; % }}}
+		elseif strcmpi(bas.name,'dronningMaud'), % {{{
+			md.basalforcings=linearbasalforcings(md);
+			md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
+			md.basalforcings.deepwater_melting_rate=3;
+			md.basalforcings.deepwater_elevation=-800;
+			md.basalforcings.upperwater_elevation=-300; % }}}
+		elseif strcmpi(bas.name,'slessor'), % {{{
+			md.basalforcings=linearbasalforcings(md);
+			md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
+			md.basalforcings.deepwater_melting_rate=1.4;
+			md.basalforcings.deepwater_elevation=-1200;
+			md.basalforcings.upperwater_elevation=-500; % }}}
+		elseif strcmpi(bas.name,'victoria'), % {{{
+			md.basalforcings=linearbasalforcings(md);
+			md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
+			md.basalforcings.deepwater_melting_rate=10;
+			md.basalforcings.deepwater_elevation=-1000;
+			md.basalforcings.upperwater_elevation=-200; % }}}
+		elseif strcmpi(bas.name,'adelie'), % {{{
+			md.basalforcings=linearbasalforcings(md);
+			md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
+			md.basalforcings.deepwater_melting_rate=4;
+			md.basalforcings.deepwater_elevation=-1000;
+			md.basalforcings.upperwater_elevation=-200; % }}}
+		elseif strcmpi(bas.name,'wilkes'), % {{{
+			md.basalforcings=linerbasalforcings(md);
+			md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
+			md.basalforcings.deepwater_melting_rate=14;
+			md.basalforcings.deepwater_elevation=-1500;
+			md.basalforcings.upperwater_elevation=-50; % }}}
+		elseif strcmpi(bas.name,'wilkins'), % {{{
+			md.basalforcings=linearbasalforcings(md);
+			md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
+			md.basalforcings.deepwater_melting_rate=6;
+			md.basalforcings.deepwater_elevation=-600;
+			md.basalforcings.upperwater_elevation=-50; % }}}
+		elseif strcmpi(bas.name,'queenmary'), % {{{
+			md.basalforcings=linearbasalforcings(md);
+			md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
+			md.basalforcings.deepwater_melting_rate=10;
+			md.basalforcings.deepwater_elevation=-1200;
+			md.basalforcings.upperwater_elevation=-200; % }}}
+		elseif strcmpi(bas.name,'mariebyrd'), % {{{
+			md.basalforcings=linearbasalforcings(md);
+			md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
+			md.basalforcings.deepwater_melting_rate=8;
+			md.basalforcings.deepwater_elevation=-700;
+			md.basalforcings.upperwater_elevation=-200; % }}}
+		elseif strcmpi(bas.name,'bellingshausen'), % {{{
+			md.basalforcings=linearbasalforcings(md);
+			md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
+			md.basalforcings.deepwater_melting_rate=7;
+			md.basalforcings.deepwater_elevation=-550;
+			md.basalforcings.upperwater_elevation=-150; % }}}
+		elseif strcmpi(bas.name,'palmer'), % {{{
+			md.basalforcings=linearbasalforcings(md);
+			md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
+			md.basalforcings.deepwater_melting_rate=.1;
+			md.basalforcings.deepwater_elevation=-1000;
+			md.basalforcings.upperwater_elevation=-50; % }}}
+		elseif strcmpi(bas.name,'peninsula'), % {{{
+			md.basalforcings=linearbasalforcings(md);
+			md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
+			md.basalforcings.deepwater_melting_rate=2;
+			md.basalforcings.deepwater_elevation=-550;
+			md.basalforcings.upperwater_elevation=-50; % }}}
+		else
+			%melt rates: {{{
+			md.basalforcings.floatingice_melting_rate=zeros(md.mesh.numberofvertices,1);
+			md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
+			load('./Data/monthly_mean_melt_adj09.mat'); fw2=fw_flux;
+			load('./Data/monthly_mean_melt_adj04.mat'); fw1=fw_flux;
+			fw_flux(:,:,end+1:end+size(fw2,3))=fw2;
+	load('./Data/grid_info_cp09.mat')
+	[xi,yi]= ll2xy(lat,lon,-1,0,71); index=delaunay(xi,yi);
+	mrate_mesh=InterpFromMeshToMesh2d(index,xi(:),yi(:),-1*reshape(squeeze(mean(fw_flux,3)),prod(size(xi)),1),md.mesh.x,md.mesh.y,'default',0);
+	pos=find(md.mask.groundedice_levelset<0 & mrate_mesh~=0); 
+	md.basalforcings.floatingice_melting_rate(pos)=mrate_mesh(pos); 
+	% }}}
+end
+md.thermal.spctemperature=[md.initialization.temperature;1]; %impose observed temperature on surface
+load(heatfluxpath)
+md.basalforcings.geothermalflux=InterpFromGridToMesh(x1,y1,bheatflx_fox,md.mesh.x,md.mesh.y,0);
+% }}}
+elseif bas.iscontinentany('greenland'),% {{{
+	%Initialize melt rates: {{{
+	md.basalforcings.floatingice_melting_rate=zeros(md.mesh.numberofvertices,1);
+	md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
+
+	%Specificy melt rates and distances:  {{{
+	%Seroussi et al., 2011 - melt rates at grounding line of 79N are about 50 m/yr
+	%Rignot and Steffen, 2008 - Petermann high rates up to 20 km out
+	if strcmpi(bas.name,'negis'), 
+		iceshelves={'gl_zac_1996','gl_79n_1996','gl_StorBistrup96'};
+		distances=[10,10,10];
+		meltrates=[7,2,4];
+		meltratesgl=[25,22,4]; %35,35
+	elseif strcmpi(bas.name,'ostenfeld'), 
+		distances=[10,10,10,10,10];
+		iceshelves={'gl_Ryder96','gl_Ostenfeld96','gl_Harder96','gl_Steensby96','gl_Brikkerne96'};
+		meltrates=[7,11,4,5,4];
+		meltratesgl=[25,26,4,5,4];
+	elseif strcmpi(bas.name,'petermann'), 
+		distances=[20,10];
+		iceshelves={'gl_Peter96','gl_Humb96'};
+		meltrates=[6,4]; %5
+		meltratesgl=[22,4]; %26
+	elseif strcmpi(bas.name,'academy'), 
+		distances=[10];
+		iceshelves={'gl_Hagen96'};
+		meltrates=[4];
+		meltratesgl=[4];
+	elseif strcmpi(bas.name,'jakobshavn'), 
+		distances=[10];
+		iceshelves={'gl_Jak'};
+		meltrates=[4];
+		meltratesgl=[109]; %Prescott et al, 2003
+	else 
+		iceshelves={};
+	end
+	%}}}
+
+	for i=1:length(iceshelves),
+		iceshelf=iceshelves{i};
+		distance=distances(i);
+		meltrate=meltrates(i);
+		meltrategl=meltratesgl(i);
+
+		in= ContourToNodes(md.mesh.x,md.mesh.y,[shppath  '/' iceshelf '_closed_polygon.shp'],1);
+		segs=MeshProfileIntersection(md.mesh.elements,md.mesh.x,md.mesh.y,[shppath  '/' iceshelf '_closed_polygon.shp']);
+		pos=find(~isnan(segs(:,1)));
+		segs=[segs(pos,1:2); segs(pos,3:4)];
+
+		%figure out nodes within ~10 km of grounding line:
+		v=find(in);
+		glpt=zeros(md.mesh.numberofvertices,1);
+		a=zeros(md.mesh.numberofvertices,1);
+		for i=v'
+			x=md.mesh.x(i); y=md.mesh.y(i);
+			dist=sqrt((segs(:,1)-x).^2+(segs(:,2)-y).^2);
+			if min(dist)<distance*1000; glpt(i)=1; else glpt(i)=0; end
+		end
+
+		md.basalforcings.floatingice_melting_rate(find(in))=meltrate;
+		md.basalforcings.floatingice_melting_rate(find(in&glpt))=meltrategl;
+
+	end
+
+	%Expand to the rest of the basin, in case of retreat of the glaciers:
+	if ~isempty(iceshelves),
+		nomelt=md.basalforcings.floatingice_melting_rate==0;
+		pos=find(nomelt);  pos2=find(~nomelt);
+		if ~isempty(pos2),
+			if length(pos2)>3,
+				md.basalforcings.floatingice_melting_rate(pos)=griddata(md.mesh.x(pos2),md.mesh.y(pos2),md.basalforcings.floatingice_melting_rate(pos2),md.mesh.x(pos),md.mesh.y(pos),'nearest');
+			end
+		end
+	end
+	%}}}
+	%Thermal spcs and heatflux:  %{{{
+	md.thermal.spctemperature=md.initialization.temperature; %impose observed temperature on surface
+
+	load(seariseg_heatflux)
+	[xi,yi]= ll2xy(md.mesh.lat,md.mesh.long,+1,39,71);
+	md.basalforcings.geothermalflux=InterpFromGridToMesh(x1,y1,bheatflx,xi,yi,0);
+	%}}}
+end % }}}
+%}}}
+%Slr: {{{
+if bas.iscontinentany('antarctica'),
+
+	delH=textread([modeldatapath '/AdhikariSciAdvTrends/AIS_delH_trend.txt']);
+	longAIS=delH(:,1);
+	latAIS=delH(:,2);
+	delHAIS=delH(:,3);
+	index=delaunay(latAIS,longAIS);
+
+	lat=md.mesh.lat; 
+	long=md.mesh.long+360; 
+	pos=find(long>360);long(pos)=long(pos)-360;
+
+	delHAIS=InterpFromMesh2d(index,longAIS,latAIS,delHAIS,long,lat);
+	northpole=find_point(md.mesh.long,md.mesh.lat,0,90); delHAIS(northpole)=0;
+
+	md.slr.deltathickness=delHAIS(md.mesh.elements)*[1;1;1]/3/100;
+else
+	delH=textread([modeldatapath '/AdhikariSciAdvTrends/GIS_delH_trend.txt']);
+	longGIS=delH(:,1);
+	latGIS=delH(:,2);
+	delHGIS=delH(:,3);
+
+	index=delaunay(latGIS,longGIS);
+
+	lat=md.mesh.lat; 
+	long=md.mesh.long+360; 
+	pos=find(long>360);long(pos)=long(pos)-360;
+
+	delHGIS=InterpFromMesh2d(index,longGIS,latGIS,delHGIS,long,lat);
+	md.slr.deltathickness=delHGIS(md.mesh.elements)*[1;1;1]/3/100;
+end	
+md.slr.sealevel=zeros(md.mesh.numberofvertices,1);
+md.slr.steric_rate=zeros(md.mesh.numberofvertices,1);
+md.slr.spcthickness=NaN*ones(md.mesh.numberofvertices,1);
+md.slr.Ngia=InterpFromMesh2d(indexGIA,longGIA,latGIA,Ngia,md.mesh.long,md.mesh.lat);
+md.slr.Ugia=InterpFromMesh2d(indexGIA,longGIA,latGIA,Ugia,md.mesh.long,md.mesh.lat);
+%}}}
+% material properties: {{{
+disp('      creating flow law paramter');
+md.materials.rheology_B=paterson(md.initialization.temperature);
+pos=find(md.materials.rheology_B<=0);
+md.materials.rheology_B(pos)=-md.materials.rheology_B(pos);
+md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
+%}}}
+%diverse: {{{
+md.miscellaneous.name=bas.name;
+% }}}
+
+	sl.icecaps{ind}=md;
+end
+%}}}
+% ParameterizeContinents {{{
+
+%glacier load: {{{ 
+delH=textread([modeldatapath '/AdhikariSciAdvTrends/GLA_delH_trend_15regions.txt']);
+longglaciers=delH(:,1); latglaciers=delH(:,2); delHglaciers=sum(delH(:,3:end),2); 
+indexglaciers=delaunay(latglaciers,longglaciers);
+%}}}
+
+for ind=sl.basinindx('continent',{'southamerica','northamerica','australia','eurasia','pacific'}),
+	disp(sprintf('Masks for basin %s\n', sl.icecaps{ind}.miscellaneous.name));
+	md=sl.icecaps{ind}; bas=sl.basins{ind}; 
+
+	%recover lat,long: 
+	[md.mesh.long,md.mesh.lat]=gdaltransform(md.mesh.x,md.mesh.y,sprintf('EPSG:%i',md.mesh.epsg),'EPSG:4326'); 
+
+	%interpolate glacier loads onto mesh:  %{{{
+	lat=md.mesh.lat; long=md.mesh.long+360; pos=find(long>360);long(pos)=long(pos)-360;
+	mdelHglaciers=InterpFromMesh2d(indexglaciers,longglaciers,latglaciers,delHglaciers,long,lat);
+	mdelHglacierse=mdelHglaciers(md.mesh.elements)*[1;1;1]/3;
+	pos=find(~md.private.bamg.landmask); mdelHglacierse(pos)=0;
+	% }}}
+	%mask:  %{{{
+	%new type of mask: 
+	md.mask=maskpsl;
+	%Figure out mask from initial mesh: deal with land and ocean masks (land include grounded ice).  %{{{
+	%first, transform land element  mask into vertex driven one.  
+	land=md.private.bamg.landmask;
+	land_mask=-ones(md.mesh.numberofvertices,1);
+
+	landels=find(land);
+	land_mask(md.mesh.elements(landels,:))=1;
+
+	%gothrough edges of each land element
+	connectedels=md.mesh.elementconnectivity(landels,:);
+	connectedisonocean=~land(connectedels);
+	sumconnectedisonocean=sum(connectedisonocean,2);
+
+	%figure out which land elements are connected to the ocean: 
+	landelsconocean=landels(find(sumconnectedisonocean));
+
+	ind1=[md.mesh.elements(landelsconocean,1); md.mesh.elements(landelsconocean,2); md.mesh.elements(landelsconocean,3)];
+	ind2=[md.mesh.elements(landelsconocean,2); md.mesh.elements(landelsconocean,3); md.mesh.elements(landelsconocean,1)];
+
+	h=waitbar(0,'Starting mask processing');
+
+	%edge ind1 and ind2: 
+	for i=1:length(ind1),
+		if ~mod(i,100),
+			waitbar(i/length(ind1),h,sprintf('%.2i%% along...',i/length(ind1)*100));
+		end
+
+		els1=md.mesh.vertexconnectivity(ind1(i),1: md.mesh.vertexconnectivity(ind1(i),end));
+		els2=md.mesh.vertexconnectivity(ind2(i),1: md.mesh.vertexconnectivity(ind2(i),end));
+		els=intersect(els1,els2);
+
+		if length(find(land(els)))==1,
+			%this edge is on the beach, 0 the edge: 
+			land_mask(ind1(i))=0;
+			land_mask(ind2(i))=0;
+		end
+	end
+	md.mask.land_levelset=land_mask;
+	close(h);
+	%}}}
+	%work on ocean, glaciers and ice:  {{{
+	%ocean: opposite of land: 
+	md.mask.ocean_levelset=-md.mask.land_levelset;
+
+	%initliaze: 
+	md.mask.ice_levelset=ones(md.mesh.numberofvertices,1); 
+
+	%%cross check that whereever we have an ice load, the mask is <0 on each vertex: if it's not th ecase, zero the ice load: 
+	pos=find(mdelHglacierse);
+	vertices=md.mesh.elements(pos,:);   %potentially having an ice load. 
+	pos=find(md.mask.land_levelset(vertices)>=0); 
+	icevertices=vertices(pos);
+	md.mask.ice_levelset(icevertices)=-md.mask.land_levelset(icevertices);
+
+	%then take care of glaciers:  don't! this is determined by the land  mask.
+	%pos=find(md.slr.deltathickness); md.mask.ice_levelset(md.mesh.elements(pos,:))=-1;
+	%pos=find(mdelHglacierse); md.mask.ice_levelset(md.mesh.elements(pos,:))=-1;
+
+
+	%grounded ice: 
+	md.mask.groundedice_levelset=-ones(md.mesh.numberofvertices,1); 
+	pos=find(mdelHglaciers); md.mask.groundedice_levelset(pos)=1;
+
+
+	%now, glaciers: 
+	md.mask.glacier_levelset=zeros(md.mesh.numberofvertices,1);
+	pos=find(mdelHglacierse); md.mask.glacier_levelset(md.mesh.elements(pos,:))=1;
+	% }}}
+	%}}}
+	%slr loading/calibration:  {{{
+	%load glaciers and ice caps from GRACE: 
+	md.slr.deltathickness=zeros(md.mesh.numberofelements,1);
+	pos=find(mdelHglacierse); md.slr.deltathickness(pos)=mdelHglacierse(pos)/100;
+
+	%initialize: 
+	md.slr.sealevel=zeros(md.mesh.numberofvertices,1);
+
+	%love numbers: 
+	md.slr.steric_rate=(1.1+.38)*ones(md.mesh.numberofvertices,1); %steric + water storage.
+	md.slr.ocean_area_scaling=1;
+	md.slr.Ngia=InterpFromMesh2d(indexGIA,longGIA,latGIA,Ngia,md.mesh.long,md.mesh.lat);
+	md.slr.Ugia=InterpFromMesh2d(indexGIA,longGIA,latGIA,Ugia,md.mesh.long,md.mesh.lat);
+	%}}}
+	%geometry:  {{{
+	di=md.materials.rho_ice/md.materials.rho_water;
+
+	md.slr.spcthickness=ones(md.mesh.numberofvertices+1,2);
+	%time tags: make sure they are far in the past and future.
+	t1=-10000; t2=+10000; md.slr.spcthickness(end,1)=-10000;
+	md.slr.spcthickness(end,2)=+10000;
+
+	%impart the glacier load: 
+	pos=find(mdelHglaciers);
+	md.slr.spcthickness(pos,2)=md.slr.spcthickness(pos,1)+(t2-t1)*mdelHglaciers(pos)/100*di;
+	%mdelHglaciers is cm/yr water equivalent
+
+	md.geometry.thickness=ones(md.mesh.numberofvertices,1);
+	md.geometry.surface=(1-di)*zeros(md.mesh.numberofvertices,1);
+	md.geometry.base=md.geometry.surface-md.geometry.thickness;
+	md.geometry.bed=md.geometry.base;
+	% }}}
+	%materials:  {{{
+	md.initialization.temperature=273.25*ones(md.mesh.numberofvertices,1);
+	md.materials.rheology_B=paterson(md.initialization.temperature);
+	md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
+	% }}}
+	%Solution: {{{
+	md.transient.issmb=0; md.transient.isstressbalance=0;
+	md.transient.ismasstransport=0; md.transient.isthermal=0;
+	md.transient.isslr=0; md.slr.loop_increment=300;
+	% }}}
+	sl.icecaps{ind}=md;
+end
+% }}}
+%Assemble Earth in 3D {{{
+
+%parameters: 
+plotting=1;
+tolerance=100;
+loneedgesdetect=0; 
+
+%create earth model by concatenating all the icecaps in 3d: 
+sl.caticecaps('tolerance',tolerance,'loneedgesdetect',loneedgesdetect);
+
+%figure out how each icecap's mesh connects to the larger earth mesh: 
+sl.intersections();
+
+%figure out connectivity: 
+disp('Mesh connectivity');
+sl.earth.mesh.vertexconnectivity=NodeConnectivity(sl.earth.mesh.elements,sl.earth.mesh.numberofvertices);
+
+%areas:
+disp('Mesh nodal areas');
+sl.earth.mesh.area=averaging(sl.earth,GetAreas3DTria(sl.earth.mesh.elements,sl.earth.mesh.x,sl.earth.mesh.y,sl.earth.mesh.z),4);
+
+%transfer a list of fields from each icecap and continent back to Earth: 
+sl.transfer('mask.groundedice_levelset');
+sl.transfer('mask.ice_levelset');
+sl.transfer('mask.ocean_levelset');
+sl.transfer('mask.land_levelset');
+sl.transfer('mask.glacier_levelset');
+sl.transfer('geometry.thickness');
+sl.transfer('geometry.surface');
+sl.transfer('geometry.base');
+sl.transfer('geometry.bed');
+sl.transfer('initialization.temperature');
+sl.transfer('materials.rheology_B');
+sl.transfer('mesh.lat');
+sl.transfer('mesh.long');
+sl.transfer('slr.deltathickness');
+sl.transfer('slr.Ngia');
+sl.transfer('slr.Ugia');
+
+%radius: 
+sl.earth.mesh.r=sqrt(sl.earth.mesh.x.^2+sl.earth.mesh.y.^2+sl.earth.mesh.z.^2);
+
+%check on the mesh transitions: {{{
+if plotting,
+	flags=ones(sl.earth.mesh.numberofvertices,1);
+	for i=1:length(sl.transitions),
+		flags(sl.transitions{i})=i;
+	end
+	plotmodel(sl.earth,'data',flags,'shading','faceted','coastline','on','coast_color','g')
+end
+%}}}}
+
+% }}}
+%Solve Sea-level equation on Earth only:  {{{
+md=sl.earth; %we don't do computations on ice sheets or land.
+
+%elastic loading from love numbers: 
+nlov=1001;
+md.slr.love_h = love_numbers('h'); md.slr.love_h(nlov+1:end)=[];
+md.slr.love_k = love_numbers('k'); md.slr.love_k(nlov+1:end)=[];
+md.slr.love_l = love_numbers('l'); md.slr.love_l(nlov+1:end)=[];
+md.slr.ocean_area_scaling = 0; 
+
+%Miscellaneous
+md.miscellaneous.name='test2004';
+
+%New stuff
+md.slr.spcthickness = NaN(md.mesh.numberofvertices,1);
+md.slr.Ngia = zeros(md.mesh.numberofvertices,1);
+md.slr.Ugia = zeros(md.mesh.numberofvertices,1);
+
+%Solution parameters
+md.slr.reltol=NaN;
+md.slr.abstol=1e-3;
+md.slr.geodetic=1;
+
+%eustatic + rigid + elastic run: 
+md.slr.rigid=1; md.slr.elastic=1; md.slr.rotation=0; 
+md.cluster=generic('name',oshostname(),'np',3);
+%md.verbose=verbose('111111111');
+md=solve(md,'Sealevelrise');
+SnoRotation=md.results.SealevelriseSolution.Sealevel;
+
+%eustatic + rigid + elastic + rotation run: 
+md.slr.rigid=1; md.slr.elastic=1; md.slr.rotation=1;
+md.cluster=generic('name',oshostname(),'np',3);
+%md.verbose=verbose('111111111');
+md=solve(md,'Sealevelrise');
+SRotation=md.results.SealevelriseSolution.Sealevel;
+
+%Fields and tolerances to track changes
+field_names     ={'noRotation','Rotation'};
+field_tolerances={1e-13,1e-13};
+field_values={SnoRotation,SRotation};
+
+%}}}
