diff --git a/examples/Helheim/Greenland.par b/examples/Helheim/Greenland.par
new file mode 100644
index 000000000..156a1226f
--- /dev/null
+++ b/examples/Helheim/Greenland.par
@@ -0,0 +1,94 @@
+%required Data
+accpath = ['/home/ModelData/Greenland/RACMO2Accumulation/SMBGreenland/smb_1961-1990.mat'];
+velpath = ['/home/ModelData/Greenland/VelMouginot/RignotGreenland2012Vel.mat'];
+
+%WARNING: we do have a levelset here, so we must set it to one even for stressbalances...
+md.transient.ismovingfront=1;
+
+disp(' Interpolating mask');
+mask = int8(interpBedmachineGreenland(md.mesh.x,md.mesh.y,'mask'));
+md.mask.ice_levelset= -1*ones(md.mesh.numberofvertices,1);
+pos = find(mask<1);
+md.mask.ice_levelset(pos)=1;
+
+disp(' reading MC bed (assumes no floating ice)');
+md.geometry.bed = interpBedmachineGreenland(md.mesh.x,md.mesh.y,'bed');
+md.geometry.base = md.geometry.bed;
+
+disp(' reading Howat surface');
+md.geometry.surface=interpBedmachineGreenland(md.mesh.x,md.mesh.y,'surface');
+pos = find(md.mask.ice_levelset>0);
+md.geometry.surface(pos) = md.geometry.base(pos)+10; %Minimum thickness
+
+md.geometry.thickness = md.geometry.surface - md.geometry.bed;
+pos=find(md.geometry.thickness<=10);
+md.geometry.surface(pos) = md.geometry.base(pos)+10; %Minimum thickness
+md.geometry.thickness = md.geometry.surface - md.geometry.bed;
+
+md.masstransport.min_thickness = 10;
+
+disp(' Adjusting ice mask');
+%Tricky part here: we want to offset the mask by one element so that we don't end up with a cliff at the transition
+pos = find(max(md.mask.ice_levelset(md.mesh.elements),[],2)>0);
+md.mask.ice_levelset(md.mesh.elements(pos,:)) = 1;
+% For the region where surface is NaN, set thickness to small value (consistency requires >0)
+pos=find((md.mask.ice_levelset<0).*(md.geometry.surface<0));
+md.mask.ice_levelset(pos)=1;
+pos=find((md.mask.ice_levelset<0).*(isnan(md.geometry.surface)));
+md.mask.ice_levelset(pos)=1;
+
+disp(' -- reconstruct thickness');
+md.geometry.thickness=md.geometry.surface-md.geometry.base;
+
+disp(' reading velocities ');
+[md.inversion.vx_obs md.inversion.vy_obs]=interpJoughinCompositeGreenland(md.mesh.x,md.mesh.y);
+pos=find(isnan(md.inversion.vx_obs) | isnan(md.inversion.vy_obs));
+md.inversion.vx_obs(pos)=0;
+md.inversion.vy_obs(pos)=0;
+md.inversion.vel_obs = sqrt(md.inversion.vx_obs.^2+md.inversion.vy_obs.^2);
+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;
+
+%drag md.drag or stress
+md.friction.coefficient = 20*ones(md.mesh.numberofvertices,1); %q = 1.
+md.friction.p = ones(md.mesh.numberofelements,1);
+md.friction.q = ones(md.mesh.numberofelements,1);
+
+%No friction on PURELY ocean element
+pos_e = find(min(md.mask.ice_levelset(md.mesh.elements),[],2)<0);
+flags=ones(md.mesh.numberofvertices,1);
+flags(md.mesh.elements(pos_e,:))=0;
+md.friction.coefficient(find(flags))=0;
+
+%flow law
+disp(' Creating flow law parameters (assume ice is at 0°C for now)');
+md.materials.rheology_n = 3*ones(md.mesh.numberofelements,1);
+md.materials.rheology_B = cuffey(273.15 - 10)*ones(md.mesh.numberofvertices,1);
+
+md.basalforcings.groundedice_melting_rate = zeros(md.mesh.numberofvertices,1);
+md.basalforcings.floatingice_melting_rate = zeros(md.mesh.numberofvertices,1);
+
+disp(' Geothermal flux from Shapiro et al.');
+md.basalforcings.geothermalflux=interpSeaRISE(md.mesh.x,md.mesh.y,'bheatflx');
+
+disp(' Setting up thermal model');
+md.initialization.temperature=min(0,interpSeaRISE(md.mesh.x,md.mesh.y,'surftemp'))+273.15;
+md.initialization.waterfraction=zeros(md.mesh.numberofvertices,1);
+md.initialization.watercolumn=zeros(md.mesh.numberofvertices,1);
+md.thermal.spctemperature=md.initialization.temperature;
+md.thermal.isenthalpy=1;
+md.thermal.isdynamicbasalspc=1;
+
+%Deal with boundary conditions:
+disp(' Set 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);
+md.stressbalance.loadingforce=0*ones(md.mesh.numberofvertices,3);
+pos=find((md.mask.ice_levelset<0).*(md.mesh.vertexonboundary));
+md.stressbalance.spcvx(pos)=md.initialization.vx(pos);
+md.stressbalance.spcvy(pos)=md.initialization.vy(pos);
+md.stressbalance.spcvz(pos)=0;
diff --git a/examples/Helheim/Helheim_tut.exp b/examples/Helheim/Helheim_tut.exp
new file mode 100644
index 000000000..20ba131cd
--- /dev/null
+++ b/examples/Helheim/Helheim_tut.exp
@@ -0,0 +1,20 @@
+## Name:./Exp/Domain_TEMP.exp
+## Icon:0
+# Points Count Value
+14 1.000000
+# X pos Y pos
+251808.9599315871 -2519299.1448521991
+232017.0235357451 -2526258.4544728124
+232815.0354161858 -2578466.8389138547
+256982.9213724715 -2601667.2788287275
+302263.2150846492 -2594953.0557058132
+309967.1091796352 -2586747.6107035289
+315046.5997917140 -2581507.9891960463
+313861.3853155622 -2571621.9108800408
+310644.3745945790 -2567964.0619031191
+309120.5274109554 -2557188.2365386733
+306496.1239280481 -2552245.1973806708
+305818.8585131043 -2537910.3838224630
+285125.3861979210 -2517350.5381584275
+251808.9599315871 -2519299.1448521991
+
diff --git a/examples/Helheim/runme.m b/examples/Helheim/runme.m
new file mode 100644
index 000000000..a28cc8b45
--- /dev/null
+++ b/examples/Helheim/runme.m
@@ -0,0 +1,86 @@
+% runme script to perform a stress balance inversion for basal friction
+% on Helheim Glacier (SE Greenland)
+
+% Step 1: Set up model domain outline and mesh
+% Step 2: Set parameters
+% Step 3: Inversion for basal friction
+
+% Note that Steps 1 and 2 require previously downloading data for velocities,
+% surface elevation and bed elevation.
+
+steps = [1:3]; % Choose which steps to run here
+
+cluster=generic('name',oshostname(),'np',8);
+
+org=organizer('repository',['./Models'],'prefix',['Model_Helheim_'],'steps',steps); clear steps;
+
+if perform(org,'Mesh'),% {{{
+
+ % Initialize uniform mesh
+ md=triangle(model,['Helheim_tut.exp'],500);
+
+ [velx vely]=interpJoughinCompositeGreenland(md.mesh.x,md.mesh.y);
+ vel = sqrt(velx.^2+vely.^2);
+
+ % Refine mesh based on surface velocities
+ md=bamg(md,'hmin',100,'hmax',1500,'field',vel,'err',5);
+ [md.mesh.lat,md.mesh.long] = xy2ll(md.mesh.x,md.mesh.y,+1,45,70);
+ md.mesh.epsg=3413;
+
+ savemodel(org,md);
+end %}}}
+if perform(org,'Param'),% {{{
+
+ md=loadmodel(org,'Mesh');
+ md=setflowequation(md,'SSA','all'); % Set flow law to SSA
+
+ md=setmask(md,'','');
+ md=parameterize(md,'Greenland.par');
+ md.miscellaneous.name = 'Helheim';
+
+ savemodel(org,md);
+end%}}}
+if perform(org,'Inversion_drag'),% {{{
+
+ md=loadmodel(org,'Param');
+
+ %Control general
+ md.inversion=m1qn3inversion(md.inversion);
+ md.inversion.iscontrol=1;
+ md.verbose=verbose('solution',false,'control',true);
+ md.transient.amr_frequency = 0;
+
+ %Cost functions
+ md.inversion.cost_functions=[101 103 501];
+ md.inversion.cost_functions_coefficients=zeros(md.mesh.numberofvertices,numel(md.inversion.cost_functions));
+ md.inversion.cost_functions_coefficients(:,1)=5000;
+ md.inversion.cost_functions_coefficients(:,2)=10;
+ md.inversion.cost_functions_coefficients(:,3)=2*50^-3;
+ pos=find(md.mask.ice_levelset>0);
+ md.inversion.cost_functions_coefficients(pos,1:2)=0;
+
+ %Controls
+ md.inversion.control_parameters={'FrictionCoefficient'};
+ md.inversion.maxsteps=100;
+ md.inversion.maxiter =100;
+ md.inversion.min_parameters=0.05*ones(md.mesh.numberofvertices,1);
+ md.inversion.max_parameters=4000*ones(md.mesh.numberofvertices,1);
+ md.inversion.control_scaling_factors=1;
+
+ %Additional parameters
+ md.stressbalance.restol=0.01;
+ md.stressbalance.reltol=0.1;
+ md.stressbalance.abstol=NaN;
+
+ %Go solve
+ md.cluster=cluster;
+
+ md=solve(md,'sb');
+
+ %Put results back into the model
+ md.friction.coefficient=md.results.StressbalanceSolution.FrictionCoefficient;
+ md.initialization.vx=md.results.StressbalanceSolution.Vx;
+ md.initialization.vy=md.results.StressbalanceSolution.Vy;
+
+ savemodel(org,md);
+end%}}}
\ No newline at end of file
diff --git a/examples/HelheimSHAKTI/runme.m b/examples/HelheimSHAKTI/runme.m
new file mode 100644
index 000000000..6c11dcaea
--- /dev/null
+++ b/examples/HelheimSHAKTI/runme.m
@@ -0,0 +1,72 @@
+% Coupled SHAKTI-ISSM simulation of Helheim, starting from previous inversion
+clear all;close all
+
+% Load Helheim Glacier model
+load Models/Model_Helheim_Inversion_drag.mat
+
+% Turn off inversion
+md.inversion.iscontrol=0;
+
+% HYDROLOGY SPECIFIC PARAMETERIZATION:
+% Change hydrology class to SHAKTI model
+md.hydrology=hydrologyshakti();
+
+% Define distributed englacial input to the subglacial system (m/yr)
+md.hydrology.englacial_input = .0*ones(md.mesh.numberofvertices,1);
+
+% Define initial water head such that water pressure is 80% of ice overburden pressure
+md.hydrology.head = 0.8*md.materials.rho_ice/md.materials.rho_freshwater*md.geometry.thickness + md.geometry.base;
+
+% Initial subglacial gap height of 0.01m everywhere
+md.hydrology.gap_height = 0.01*ones(md.mesh.numberofelements,1);
+
+% Typical bed bump bump spacing -- must be non-zero
+md.hydrology.bump_spacing = 1.0*ones(md.mesh.numberofelements,1);
+
+% Typical bed bump height -- set to zero to turn off opening by sliding
+md.hydrology.bump_height = 0.0*ones(md.mesh.numberofelements,1);
+
+% Initial Reynolds number (start at Re=1000 everywhere)
+md.hydrology.reynolds= 1000*ones(md.mesh.numberofelements,1);
+
+% Boundary conditions
+md.hydrology.spchead = NaN(md.mesh.numberofvertices,1);
+
+% Set head=0 for thin ice, everywhere <=10m ice thickness
+pos=find(md.geometry.thickness==10);
+md.hydrology.spchead(pos)=0;
+
+% Set pressure in fjord equal to hydrostatic pressure of fjord water
+pos=find(md.mask.ice_levelset>0);
+md.hydrology.spchead(pos)=0;
+
+% % Deal with terminus b.c. in SHAKTI
+% terminus=ContourToNodes(md.mesh.x,md.mesh.y,'Exp/terminus_big.exp',1);
+% pos=find(terminus & md.mask.ice_levelset<=0);
+% md.hydrology.spchead(pos)=0;
+
+md.hydrology.moulin_input = zeros(md.mesh.numberofvertices,1); % No moulin inputs
+md.hydrology.neumannflux=zeros(md.mesh.numberofelements,1); % No-flux b.c. on domain boundary
+
+% Coupling and friction
+md.transient=deactivateall(md.transient);
+md.transient.isstressbalance=1; % 1=Solve for ice velocity, 0=turn off stress balance
+md.transient.ishydrology=1;
+
+% Friction
+Neff = md.materials.rho_ice*md.constants.g*md.geometry.thickness-md.materials.rho_water*md.constants.g*(md.hydrology.head - md.geometry.base); %Initial effective pressure
+md.friction.effective_pressure=Neff;
+md.friction.coupling = 4; % Change this to couple (4) /uncouple with prescribed N (3) ***
+
+% Specify that you want to run the model on your current computer
+% Change the number of processors according to your machine
+md.cluster=generic('np',8);
+
+% Define the time stepping scheme
+md.timestepping.time_step=1*3600/md.constants.yts; % Time step (in years)
+md.timestepping.final_time=30/365; % Final time (in years)
+md.settings.output_frequency=24; % To save model output less frequently than every time step
+
+md.verbose.solution=1;
+
+md=solve(md,'Transient');