Index: issm/trunk-jpl/examples/MeltingSensitivity/DomainOutline.exp
===================================================================
--- issm/trunk-jpl/examples/MeltingSensitivity/DomainOutline.exp	(revision 20735)
+++ 	(revision )
@@ -1,101 +1,0 @@
-## Name:
-## Icon:0
-# Points Count Value
-95 25000.000000
-# X pos Y pos
--1712113.0179281300 -349656.0205056490
--1711509.1917759699 -349354.1074295690
--1709240.6586108401 -349354.1074295690
--1707080.1508345299 -348922.0058743060
--1705459.7700022999 -347193.5996532570
--1703083.2114483600 -345249.1426545770
--1701138.7544496800 -344168.8887664210
--1698438.1197292900 -343952.8379887900
--1698006.0181740201 -341900.3556012940
--1693793.0280102200 -340279.9747690600
--1691092.3932898301 -339631.8224361670
--1688391.7585694401 -339415.7716585360
--1686447.3015707601 -340063.9239914290
--1683530.6160727399 -341036.1524907690
--1682558.3875734000 -343304.6856558960
--1682134.7209454600 -344182.1487658350
--1682126.2860181299 -345681.2442098390
--1681072.0444963200 -345932.4393879640
--1680634.4718407800 -346682.5639403050
--1680571.9614614199 -348557.8753211580
--1677446.4424933300 -348370.3441830730
--1675508.6207331200 -349495.5310115840
--1676616.9911885399 -353567.0975933760
--1672991.8601872099 -353799.6175753630
--1668839.1631938200 -353999.1991486380
--1667132.2298986400 -355121.4651541420
--1666069.5534494901 -357496.8595698890
--1662944.0344814099 -358309.4945015920
--1660756.1712037399 -358747.0671571240
--1657922.9836098100 -361682.8074238540
--1656579.2228470000 -364003.8487410000
--1653004.8841628900 -364748.0635758530
--1653442.4568184200 -361622.5446077650
--1651129.5727820301 -355684.0585683980
--1648929.2928810001 -350845.9692010000
--1646260.7698359799 -347802.2804785470
--1644571.8968503401 -341785.6704672340
--1643503.3064999001 -340869.0986596620
--1643190.7546030900 -336680.9032424240
--1642690.6715682000 -335243.1645171040
--1641913.0916588283 -332341.0814454452
--1640683.9739821283 -330436.6112348017
--1637816.0327364956 -328532.1410241582
--1635241.2508976001 -328652.8685499910
--1631765.4107996500 -327017.1790921330
--1625934.5618617306 -322342.6128395668
--1620608.3852626982 -320438.1426289233
--1616101.6204481323 -320438.1426289233
--1611594.8556335662 -320914.2601815842
--1607907.5026034669 -321866.4952869060
--1602991.0318966676 -324723.2006028712
--1599713.3847588014 -327579.9059188365
--1593977.5022675355 -328532.1410241582
--1590699.8551296694 -334721.6692087495
--1589259.8391499999 -338300.0840570000
--1588627.0898328500 -343572.9950327750
--1588056.6000060199 -345432.3563530700
--1586993.9235568701 -347432.6884926460
--1589259.8391499999 -349315.9832070000
--1585493.6744521901 -352996.1122558430
--1585118.6121760199 -354496.3613605250
--1585861.9157163899 -356240.3342513760
--1586368.8197632500 -358747.0671571240
--1587889.0972700799 -360996.6006787950
--1586681.3716600600 -363935.4286441500
--1585806.2263489999 -365748.2296456410
--1586743.8820394201 -368373.6655788350
--1586056.2678664399 -370123.9562009640
--1585493.6744521901 -371124.1222707520
--1586243.7990045301 -373562.0270658610
--1414257.6683934701 -327460.5482541260
--1412220.4351378100 -325972.8942023180
--1404515.9240538401 -296165.1737279960
--1387984.9332285200 -267002.0106682320
--1360067.9373810000 -107884.1934940000
--1423274.3497420000 -81217.9259650000
--1402446.7987670000 -45916.9921070000
--1478343.8065599999 60338.8188030000
--1503760.4789370000 36334.1837800000
--1528509.5619775900 63667.6809196406
--1623055.4918765700 62841.1506965373
--1725734.5897359999 -61372.6193030000
--1719614.6457630000 -69022.5492690000
--1759088.2843859999 -114004.1374660000
--1755416.3180020000 -117064.1094530000
--1761230.2647760001 -126856.0198080000
--1758170.2927900001 -130527.9861920000
--1762148.2563720001 -135423.9413700000
--1764902.2311590000 -131445.9777880000
--1769798.1863370000 -134199.9525750000
--1775000.1387139999 -131445.9777880000
--1794277.9622269999 -153477.7760890000
--1712576.7101950001 -225081.1205660000
--1729939.5891459100 -335652.4567604720
--1712113.0179281300 -349656.0205056490
-
Index: issm/trunk-jpl/examples/MeltingSensitivity/Pig.par
===================================================================
--- issm/trunk-jpl/examples/MeltingSensitivity/Pig.par	(revision 20735)
+++ 	(revision )
@@ -1,129 +1,0 @@
-% Parameters to change/Try
-friction_coefficient = 10; % default [10]
-Temp_change          =  0;  % default [0 K]
-
-%Name and Coordinate system
-md.miscellaneous.name='PIG';
-md.mesh.epsg=3031;
-
-%NetCdf Loading
-disp('   Loading SeaRISE data from NetCDF');
-ncdata='../Data/Antarctica_5km_withshelves_v0.75.nc';
-x1    = ncread(ncdata,'x1');
-y1    = ncread(ncdata,'y1');
-usrf  = ncread(ncdata,'usrf')';
-topg  = ncread(ncdata,'topg')';
-temp  = ncread(ncdata,'presartm')';
-smb   = ncread(ncdata,'presprcp')';
-gflux = ncread(ncdata,'bheatflx_fox')';
-
-disp('   Loading velocities data from NetCDF');
-nsidc_vel='../Data/Antarctica_ice_velocity.nc';
-
-xmin = ncreadatt(nsidc_vel,'/','xmin');
-xmin = strtrim(xmin);  %this is a string, and we need to recover the double value
-xmin = xmin(1:end-2);  %get rid of the unit
-xmin = str2num(xmin);  %convert to double
-
-ymax = ncreadatt(nsidc_vel,'/','ymax'); 
-ymax = strtrim(ymax);  
-ymax = ymax(1:end-2);  
-ymax = str2num(ymax); 
-
-nx = ncreadatt(nsidc_vel,'/','nx');
-ny = ncreadatt(nsidc_vel,'/','ny');
-
-spacing = ncreadatt(nsidc_vel,'/','spacing'); 
-spacing = strtrim(spacing);  
-spacing = spacing(1:end-2);  
-spacing = str2num(spacing); 
-
-velx = double(ncread(nsidc_vel,'vx'));
-vely = double(ncread(nsidc_vel,'vy'));
-
-x2=xmin+(0:1:nx)'*spacing; 
-x2=double(x2);
-
-y2=(ymax-ny*spacing)+(0:1:ny)'*spacing; 
-y2=double(y2);
-
-%Geometry
-disp('   Interpolating surface and ice base');
-md.geometry.base    = InterpFromGridToMesh(x1,y1,topg,md.mesh.x,md.mesh.y,0);
-md.geometry.surface = InterpFromGridToMesh(x1,y1,usrf,md.mesh.x,md.mesh.y,0);
-clear usrf, topg;
-
-disp('   Constructing thickness');
-md.geometry.thickness=md.geometry.surface-md.geometry.base;
-
-%ensure hydrostatic equilibrium on ice shelf: 
-di=md.materials.rho_ice/md.materials.rho_water;
-
-%Get the node numbers of floating nodes
-pos=find(md.mask.groundedice_levelset<0); 
-
-%apply a flotation criterion on the precedingly defined nodes and
-%redefine base and thickness accordingly
-md.geometry.thickness(pos)=1/(1-di)*md.geometry.surface(pos);
-md.geometry.base(pos)=md.geometry.surface(pos)-md.geometry.thickness(pos);
-md.geometry.hydrostatic_ratio=ones(md.mesh.numberofvertices,1);
-
-%Set min thickness to 1 meter
-pos0=find(md.geometry.thickness<=0);
-md.geometry.thickness(pos0)=1;
-md.geometry.surface=md.geometry.thickness+md.geometry.base;
-
-%Initialization parameters
-disp('   Interpolating temperatures');
-md.initialization.temperature=InterpFromGridToMesh(x1,y1,temp,md.mesh.x,md.mesh.y,0)+273.15+Temp_change;
-clear temp;
-
-disp('   Set observed velocities')
-vx_obs=InterpFromGridToMesh(x2,y2,flipud(velx'),md.mesh.x,md.mesh.y,0);
-vy_obs=InterpFromGridToMesh(x2,y2,flipud(vely'),md.mesh.x,md.mesh.y,0);
-clear velx vely;
-
-vel_obs=sqrt(vx_obs.^2+vy_obs.^2);
-md.initialization.vx=vx_obs;
-md.initialization.vy=vy_obs;
-md.initialization.vz=zeros(md.mesh.numberofvertices,1);
-md.initialization.vel=vel_obs;
-
-disp('   Set Pressure');
-md.initialization.pressure=md.materials.rho_ice*md.constants.g*md.geometry.thickness;
-
-disp('   Construct ice rheological properties');
-md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
-md.materials.rheology_B=paterson(md.initialization.temperature);
-
-%Forcings
-disp('   Interpolating surface mass balance');
-md.smb.mass_balance=InterpFromGridToMesh(x1,y1,smb,md.mesh.x,md.mesh.y,0);
-md.smb.mass_balance=md.smb.mass_balance*md.materials.rho_water/md.materials.rho_ice;
-clear smb;
-
-disp('   Set geothermal heat flux');
-md.basalforcings.geothermalflux=InterpFromGridToMesh(x1,y1,gflux,md.mesh.x,md.mesh.y,0);
-clear gflux;
-
-%Friction and inversion set up
-disp('   Construct basal friction parameters');
-md.friction.coefficient=friction_coefficient*ones(md.mesh.numberofvertices,1);
-md.friction.p=ones(md.mesh.numberofelements,1);
-md.friction.q=ones(md.mesh.numberofelements,1);
-
-%no friction applied on floating ice
-pos=find(md.mask.groundedice_levelset<0);
-md.friction.coefficient(pos)=0;
-
-md.inversion=m1qn3inversion();
-md.inversion.vx_obs=vx_obs;
-md.inversion.vy_obs=vy_obs;
-md.inversion.vel_obs=vel_obs;
-
-disp('   Set boundary conditions');
-md=SetMarineIceSheetBC(md);
-md.basalforcings.floatingice_melting_rate = zeros(md.mesh.numberofvertices,1);
-md.basalforcings.groundedice_melting_rate = zeros(md.mesh.numberofvertices,1);
-md.thermal.spctemperature     = [md.initialization.temperature;1]; %impose observed temperature on surface
-md.masstransport.spcthickness    = NaN*ones(md.mesh.numberofvertices,1);
Index: issm/trunk-jpl/examples/MeltingSensitivity/runme.m
===================================================================
--- issm/trunk-jpl/examples/MeltingSensitivity/runme.m	(revision 20735)
+++ issm/trunk-jpl/examples/MeltingSensitivity/runme.m	(revision 20737)
@@ -1,221 +1,73 @@
-%Which steps to be performed
-steps=[1] ;
+step=[1];
 
-%Run Steps
+if step==1 %Transient Run #1
 
-%Mesh Generation #1
-if any(steps==1)
+	md = loadmodel('./Models/PIG.Control_drag');	
 
-	md.miscellaneous.name='PIG.Mesh_generation';
+	md.inversion.iscontrol=0;
+	md.transient.ismasstransport=1;
+	md.transient.isstressbalance=1;
+	md.transient.isgroundingline=1;
+	md.transient.ismovingfront=0;
+	md.transient.isthermal=0;
+	
+	md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
+	md.basalforcings.floatingice_melting_rate=30*ones(md.mesh.numberofvertices,1);
+	
+	md.timestepping.time_step=0.1;
+	md.timestepping.final_time=10;
+	md.transient.requested_outputs={'default','IceVolume','IceVolumeAboveFloatation'}
 
-	%Mesh parameters
-	domain =['./DomainOutline.exp'];
-	hmax=40000;    % maximum element size of the final mesh
-	hmin=5000;     % minimum element size of the final mesh
-	hinit=10000;   % element size for the initial mesh
-	gradation=1.7; % maximum size ratio between two neighboring elements
-	err=8;         % maximum error between interpolated and control field
+	md=solve(md,TransientSolutionEnum);
 
-	% Generate an initial uniform mesh (resolution = hinit m)
-	md=bamg(model,'domain',domain,'hmax',hinit,'MaxCornerAngle',1);
-
-	%ploting
-	plotmodel(md,'data','mesh')
-
-	% Load Velocities
-	% http://nsidc.org/data/nsidc-0484.html
-	nsidc_vel='../Data/Antarctica_ice_velocity.nc'; 	
-
-	% Get necessary data to build up the velocity grid
-	xmin = ncreadatt(nsidc_vel,'/','xmin');
-	xmin = strtrim(xmin);  % this is a string, and we need to recover the double value
-	xmin = xmin(1:end-2);  % get rid of the unit
-	xmin = str2num(xmin);  % convert to double
-	
-	ymax = ncreadatt(nsidc_vel,'/','ymax'); 
-	ymax = strtrim(ymax);  
-	ymax = ymax(1:end-2);  
-	ymax = str2num(ymax); 
-	
-	nx = ncreadatt(nsidc_vel,'/','nx');
-	ny = ncreadatt(nsidc_vel,'/','ny');
-	
-	spacing = ncreadatt(nsidc_vel,'/','spacing'); 
-	spacing = strtrim(spacing);
-	spacing = spacing(1:end-2);  
-	spacing = str2num(spacing); 
-	
-	% Get velocities (Note: You can use ncdisp('file') to see an ncdump)
-	vx = double(ncread(nsidc_vel,'vx'));
-	vy = double(ncread(nsidc_vel,'vy'));
-	
-	x=xmin+(0:1:nx)'*spacing; 
-	x=double(x);
-	y=(ymax-ny*spacing)+(0:1:ny)'*spacing; 
-	y=double(y);
-
-	% Interpolate velocities onto coarse mesh
-	vx_obs=InterpFromGridToMesh(x,y,flipud(vx'),md.mesh.x,md.mesh.y,0);
-	vy_obs=InterpFromGridToMesh(x,y,flipud(vy'),md.mesh.x,md.mesh.y,0);
-	vel_obs=sqrt(vx_obs.^2+vy_obs.^2);
-	clear vx vy x y;
-
-	% Adapt the mesh to minimize error in velocity interpolation
-	md=bamg(md,'hmax',hmax,'hmin',hmin,'gradation',gradation,'field',vel_obs,'err',err);
-	
-	%ploting
-	plotmodel(md,'data','mesh')
-
-	% Convert x,y coordinates (Polar stereo) to lat/lon
-	[md.mesh.lat,md.mesh.long]=xy2ll(md.mesh.x,md.mesh.y,-1);
-	
 	% Save model
-	save ./Models/PIG.Mesh_generation md;
+	save ./Models/PIG.Transient md;
 end
 
-%Masks #2
-if any(steps==2) 
+if step==2 %High Melt #2
+	md = loadmodel('./Models/PIG.Transient');	
 
-	md = loadmodel('./Models/PIG.Mesh_generation');	
+	md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
+	md.basalforcings.floatingice_melting_rate=60*ones(md.mesh.numberofvertices,1);
+	
+	md.timestepping.time_step=0.1;
+	md.timestepping.final_time=10;
+	md.transient.requested_outputs={'default','IceVolume','IceVolumeAboveFloatation'}
 
-	% Load SeaRISe dataset for Antarctica  
-	% http://websrv.cs.umt.edu/isis/index.php/Present_Day_Antarctica
-	searise='../Data/Antarctica_5km_withshelves_v0.75.nc';
-	
-	%read thickness mask from SeaRISE
-	x1=double(ncread(searise,'x1'));
-	y1=double(ncread(searise,'y1'));
-	thkmask=double(ncread(searise,'thkmask'));
-	
-	%interpolate onto our mesh vertices
-	groundedice=double(InterpFromGridToMesh(x1,y1,thkmask',md.mesh.x,md.mesh.y,0));
-	groundedice(groundedice<=0)=-1;
-	clear thkmask;
+	md=solve(md,TransientSolutionEnum);
 
-	%fill in the md.mask structure
-	md.mask.groundedice_levelset=groundedice; %ice is grounded for mask equal one
-	md.mask.ice_levelset=-1*ones(md.mesh.numberofvertices,1);%ice is present when negatvie
-
-	%ploting
-	plotmodel(md,'data',md.mask.groundedice_levelset,'title','grounded/floating','data',md.mask.ice_levelset,'title','ice/no-ice')
-	
-	% Save model
-	save ./Models/PIG.SetMask md;
+	save ./Models/PIG.HighMelt md;
 end
 
-%Parameterization #3
-if any(steps==3) 
+if step==3 %High surface mass balance #3
+	md = loadmodel('./Models/PIG.Transient');	
 
-	md = loadmodel('./Models/PIG.SetMask');
-	md = parameterize(md,'./Pig.par');
+	md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
+	md.basalforcings.floatingice_melting_rate=30*ones(md.mesh.numberofvertices,1);
 
-	% Use a MacAyeal flow model
-	md = setflowequation(md,'SSA','all');
+	md.smb.mass_balance=2*md.smb.mass_balance;
 	
-	% Save model
-	save ./Models/PIG.Parameterization md;
+	md.timestepping.time_step=0.1;
+	md.timestepping.final_time=10;
+
+	md.transient.requested_outputs={'default','IceVolume','IceVolumeAboveFloatation'}
+
+	md=solve(md,TransientSolutionEnum);
+
+	save ./Models/PIG.HighMelt md;
 end
 
-%Control Method #4
-if any(steps==4)
+if step==4 %Ice Front retreat
+	md = loadmodel('./Models/PIG.Transient');	
 
-	md = loadmodel('./Models/PIG.Parameterization');
+	md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
+	md.basalforcings.floatingice_melting_rate=30*ones(md.mesh.numberofvertices,1);
 
-	% Control general
-	md.inversion.iscontrol=1;
-	md.inversion.maxsteps=20;
-	md.inversion.maxiter=40;
-	md.inversion.dxmin=0.1;
-	md.inversion.gttol=1.0e-4;
-	md.verbose=verbose('solution',true,'control',true);
+	md.timestepping.time_step=0.1;
+	md.timestepping.final_time=10;
+	md.transient.requested_outputs={'default','IceVolume','IceVolumeAboveFloatation'}
 
-	% Cost functions
-	md.inversion.cost_functions=[101 103 501];
-	md.inversion.cost_functions_coefficients=ones(md.mesh.numberofvertices,3);
-	md.inversion.cost_functions_coefficients(:,1)=1;
-	md.inversion.cost_functions_coefficients(:,2)=1;
-	md.inversion.cost_functions_coefficients(:,3)=8e-15;
+	md=solve(md,TransientSolutionEnum);
 
-	% Controls
-	md.inversion.control_parameters={'FrictionCoefficient'};
-	md.inversion.min_parameters=1*ones(md.mesh.numberofvertices,1);
-	md.inversion.max_parameters=200*ones(md.mesh.numberofvertices,1);
-
-	% Additional parameters
-	md.stressbalance.restol=0.01;
-	md.stressbalance.reltol=0.1;
-	md.stressbalance.abstol=NaN;
-
-	% Solve
-	md.toolkits=toolkits;
-	md.cluster=generic('name',oshostname,'np',2);
-	md=solve(md,StressbalanceSolutionEnum);
-
-	% Update model friction fields accordingly
-	md.friction.coefficient=md.results.StressbalanceSolution.FrictionCoefficient;
-
-	plotmodel(md,'data',md.friction.coefficient)
-
-	% Save model
-	save ./Models/PIG.Control_drag md;
+	save ./Models/PIG.HighMelt md;
 end
-
-%Plot #5
-if any(steps==5)
-
-	md = loadmodel('./Models/PIG.Control_drag');
-
-	plotmodel(md,'nlines',2,'ncols',2,'unit#all','km','axis#all','equal',...
-		'xlim#all',[min(md.mesh.x) max(md.mesh.x)]/10^3,...
-		'ylim#all',[min(md.mesh.y) max(md.mesh.y)]/10^3,...
-		'FontSize#all',12,...
-		'data',md.initialization.vel,'title','Observed velocity',...
-		'data',md.results.StressbalanceSolution.Vel,'title','Modeled Velocity',...
-		'data',md.geometry.base,'title','Bed elevation',...
-		'data',md.results.StressbalanceSolution.FrictionCoefficient,'title','Friction Coefficient',...
-		'colorbar#all','on','colorbartitle#1-2','[m/yr]',...
-		'caxis#1-2',([1.5,4000]),...
-		'colorbartitle#3','[m]', 'log#1-2',10);
-end
-
-%HO #6
-if any(steps==6)
-
-	% Load Model
-
-	% Disable inversion
-
-	% Extrude Mesh
-
-	% Set Flowequation
-
-	% Solve
-
-	% Save Model
-
-end
-
-%Plot #7
-if any(steps==7)
-
-	mdHO = loadmodel('./Models/PIG.ModelHO');
-	mdSSA = loadmodel('./Models/PIG.Control_drag');
-
-	basal=find(mdHO.mesh.vertexonbase);
-	surf=find(mdHO.mesh.vertexonsurface);
-
-	plotmodel(mdHO,'nlines',3,'ncols',2,'unit#all','km','axis#all','equal',...
-						'xlim#all',[min(mdHO.mesh.x) max(mdHO.mesh.x)]/10^3,...
-						'ylim#all',[min(mdHO.mesh.y) max(mdHO.mesh.y)]/10^3,...
-						'FontSize#all',12,...
-						'data',mdHO.initialization.vel,'title','Observed velocity',...
-						'data',(mdHO.results.StressbalanceSolution.Vel(surf)-mdHO.initialization.vel(surf)),'title','(HO-observed) velocities',...
-						'data',mdSSA.results.StressbalanceSolution.Vel,'title','Modeled SSA Velocity',...
-						'data',(mdHO.results.StressbalanceSolution.Vel(surf)-mdSSA.results.StressbalanceSolution.Vel),'title','(HO-SSA) velocities',...
-						'data',mdHO.results.StressbalanceSolution.Vel,'title','Modeled HO surface Velocities',...
-						'data',(mdHO.results.StressbalanceSolution.Vel(surf)-mdHO.results.StressbalanceSolution.Vel(basal)),'title','(HOsurf-HO base) velocities',...
-						'caxis#1',([1.5,4000]),'caxis#3',([1.5,4000]),'caxis#5',([1.5,4000]),...
-						'colorbar#all','on','view#all',2,...
-						'colorbartitle#all','[m/yr]',...
-						'layer#5',1, 'log#1', 10,'log#3', 10,'log#5', 10);
-end
