Index: /issm/trunk-jpl/test/NightlyRun/test2051.m
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test2051.m	(revision 14849)
+++ /issm/trunk-jpl/test/NightlyRun/test2051.m	(revision 14850)
@@ -1,10 +1,30 @@
-%GIA test, inspired on test101
-md=triangle(model(),'../Exp/Square.exp',100000.);
+% Benchmark experiments (Figure A2a Ivins and James, 1999, Geophys. J. Int.) 
+md=triangle(model(),'../Exp/RoundFrontEISMINT.exp',100000.);
 md=setmask(md,'','');
-md=parameterize(md,'../Par/SquareSheetConstrained.par');
-md.timestepping.start_time=2400000; %2,400 kyr
-md.timestepping.final_time=2500000; %2,400 kyr
-md.geometry.thickness=[zeros(md.mesh.numberofvertices+1,1), [md.geometry.thickness; 1e-4], [md.geometry.thickness; md.timestepping.start_time]];
-md.cluster=generic('name',oshostname(),'np',3);
+md=parameterize(md,'../Par/GiaBenchmarksAB.par');
+
+%% indicate what you want to compute 
+md.gia.output_rates=1;           % just want "w" solution 
+md.gia.cross_section_shape=1;    % for square-edged x-section 
+
+%% define loading history 
+md.timestepping.start_time=2002100; % after 2 kyr of deglaciation 
+md.timestepping.final_time=2500000; % 2,500 kyr
+% In order to run Ivins1999 benchmarks, for now, we do things manually as follows: 
+% 1. numtimes is fixed to be five 
+% 2. final_time given above is hard-coded in "GiaDeflectionCorex.cpp"
+% 3. evaluation time = start_time 
+% 4. define new load times, according to the chosen evaluation time (care that general load history is not altered). This is really important, because of #1 and also becasue of the fact that futre loading times are not recognized.
+% 5. Ice thickness assoiated with evaluation time is actually the same at final_time
+md.geometry.thickness=[...
+	[md.geometry.thickness*0.0; 0.0],...
+	[md.geometry.thickness; 1000],...
+	[md.geometry.thickness; 2000000],...
+	[md.geometry.thickness*0.0; 2000100],...
+	[md.geometry.thickness*0.0; md.timestepping.start_time],...
+	];
+
+%% solve for GIA deflection 
+md.cluster=generic('name',oshostname(),'np',8);
 md.verbose=verbose('1111111');
 md=solve(md,GiaSolutionEnum());
