Index: /issm/trunk-jpl/test/NightlyRun/test2051.m
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test2051.m	(revision 25289)
+++ /issm/trunk-jpl/test/NightlyRun/test2051.m	(revision 25290)
@@ -1,3 +1,3 @@
-%Test Name: GiaIvinsBenchmarksAB2dA1
+%Test Name: GiaIvinsBenchmarksAB2dA
 % Benchmark experiments (Figure A2a Ivins and James, 1999, Geophys. J. Int.) 
 md=triangle(model(),'../Exp/RoundFrontEISMINT.exp',200000.);
@@ -8,7 +8,9 @@
 md.gia.cross_section_shape=1;    % for square-edged x-section 
 
-%% define loading history 
+% evaluation time (termed start_time) 
 md.timestepping.start_time=2002100; % after 2 kyr of deglaciation 
 md.timestepping.final_time=2500000; % 2,500 kyr
+
+%% define loading history {{{ 
 md.geometry.thickness=[...
 	[md.geometry.thickness*0.0; 0.0],...
@@ -18,15 +20,45 @@
 	[md.geometry.thickness*0.0; md.timestepping.start_time],...
 	];
+% }}} 
+
+% find out elements that have zero loads throughout the loading history.
+pos = find(sum(abs(md.geometry.thickness(1:end-1,:)),2)==0);
+md.mask.ice_levelset(pos)=1; % no-ice.
+
+md.cluster=generic('name',oshostname(),'np',3);
+md.verbose=verbose('1111111');
 
 %% solve for GIA deflection 
-md.cluster=generic('name',oshostname(),'np',3);
-md.verbose=verbose('1111111');
 md=solve(md,'Gia');
 
+%Test Name: GiaIvinsBenchmarksAB2dA1
+U_AB2dA1 = md.results.GiaSolution.UGia; 
+URate_AB2dA1 = md.results.GiaSolution.UGiaRate; 
+
+%Test Name: GiaIvinsBenchmarksAB2dA2
+%% different evaluation time. {{{ 
+md.timestepping.start_time=2005100; % after 5 kyr of deglaciation 
+md.geometry.thickness(end,end) = md.timestepping.start_time; 
+
+md=solve(md,'Gia');
+
+U_AB2dA2 = md.results.GiaSolution.UGia; 
+URate_AB2dA2 = md.results.GiaSolution.UGiaRate; 
+% }}} 
+
+%Test Name: GiaIvinsBenchmarksAB2dA3
+%% different evaluation time. {{{ 
+md.timestepping.start_time=2010100; % after 10 kyr of deglaciation 
+md.geometry.thickness(end,end) = md.timestepping.start_time; 
+
+md=solve(md,'Gia');
+
+U_AB2dA3 = md.results.GiaSolution.UGia; 
+URate_AB2dA3 = md.results.GiaSolution.UGiaRate; 
+% }}} 
+
 %Fields and tolerances to track changes
-field_names     ={'UGia','UGiaRate'};
-field_tolerances={1e-13,1e-13};
-field_values={...
-	(md.results.GiaSolution.UGia),...
-	(md.results.GiaSolution.UGiaRate),...
-	};
+field_names     ={'U_AB2dA1','URate_AB2dA1','U_AB2dA2','URate_AB2dA2','U_AB2dA3','URate_AB2dA3'};
+field_tolerances={1e-13,1e-13,1e-13,1e-13,1e-13,1e-13};
+field_values={U_AB2dA1,URate_AB2dA1,U_AB2dA2,URate_AB2dA2,U_AB2dA3,URate_AB2dA3}; 
+
Index: /issm/trunk-jpl/test/NightlyRun/test2052.m
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test2052.m	(revision 25289)
+++ /issm/trunk-jpl/test/NightlyRun/test2052.m	(revision 25290)
@@ -1,32 +1,63 @@
-%Test Name: GiaIvinsBenchmarksAB2dA2
+%Test Name: GiaIvinsBenchmarksAB2dC
 % Benchmark experiments (Figure A2a Ivins and James, 1999, Geophys. J. Int.) 
 md=triangle(model(),'../Exp/RoundFrontEISMINT.exp',200000.);
 md=setmask(md,'','');
-md=parameterize(md,'../Par/GiaIvinsBenchmarksAB.par');
+md=parameterize(md,'../Par/GiaIvinsBenchmarksCD.par');
 
 %% indicate what you want to compute 
 md.gia.cross_section_shape=1;    % for square-edged x-section 
 
-%% define loading history 
-md.timestepping.start_time=2005100; % after 5 kyr of deglaciation 
+% evaluation time (termed start_time) 
+md.timestepping.start_time=0.3;     % for t \approx 0 kyr : to get eleastic response!   
 md.timestepping.final_time=2500000; % 2,500 kyr
+
+%% define loading history {{{ 
 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],...
+	[md.geometry.thickness/2.0; 0.1],...
+	[md.geometry.thickness; 0.2],...
+	[md.geometry.thickness; md.timestepping.start_time],...
 	];
+% }}} 
+
+% find out elements that have zero loads throughout the loading history.
+pos = find(sum(abs(md.geometry.thickness(1:end-1,:)),2)==0);
+md.mask.ice_levelset(pos)=1; % no-ice.
+
+md.cluster=generic('name',oshostname(),'np',3);
+md.verbose=verbose('1111111');
 
 %% solve for GIA deflection 
-md.cluster=generic('name',oshostname(),'np',3);
-md.verbose=verbose('1111111');
 md=solve(md,'Gia');
 
+%Test Name: GiaIvinsBenchmarksAB2dC1
+U_AB2dC1 = md.results.GiaSolution.UGia; 
+URate_AB2dC1 = md.results.GiaSolution.UGiaRate; 
+
+%Test Name: GiaIvinsBenchmarksAB2dC2
+%% different evaluation time. {{{ 
+md.timestepping.start_time=1000.3;  % for t \approx 1 kyr 
+md.geometry.thickness(end,end) = md.timestepping.start_time; 
+
+md=solve(md,'Gia');
+
+U_AB2dC2 = md.results.GiaSolution.UGia; 
+URate_AB2dC2 = md.results.GiaSolution.UGiaRate; 
+% }}} 
+
+%Test Name: GiaIvinsBenchmarksAB2dC3
+%% different evaluation time. {{{ 
+md.timestepping.start_time=2400000; % for t \approx \infty 
+md.geometry.thickness(end,end) = md.timestepping.start_time; 
+
+md=solve(md,'Gia');
+
+U_AB2dC3 = md.results.GiaSolution.UGia; 
+URate_AB2dC3 = md.results.GiaSolution.UGiaRate; 
+% }}} 
+
 %Fields and tolerances to track changes
-field_names     ={'UGia','UGiaRate'};
-field_tolerances={1e-13,1e-13};
-field_values={...
-	(md.results.GiaSolution.UGia),...
-	(md.results.GiaSolution.UGiaRate),...
-	};
+field_names     ={'U_AB2dC1','URate_AB2dC1','U_AB2dC2','URate_AB2dC2','U_AB2dC3','URate_AB2dC3'};
+field_tolerances={1e-13,1e-13,1e-13,1e-13,1e-13,1e-13};
+field_values={U_AB2dC1,URate_AB2dC1,U_AB2dC2,URate_AB2dC2,U_AB2dC3,URate_AB2dC3}; 
+
Index: /issm/trunk-jpl/test/NightlyRun/test2052.py
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test2052.py	(revision 25289)
+++ /issm/trunk-jpl/test/NightlyRun/test2052.py	(revision 25290)
@@ -1,3 +1,3 @@
-#Test Name: GiaIvinsBenchmarksAB2dA2
+#Test Name: GiaIvinsBenchmarksAB2dC1
 from socket import gethostname
 
@@ -11,20 +11,19 @@
 
 
-#Benchmark experiments (Figure A2a Ivins and James, 1999, Geophys. J. Int.)
+#Benchmark experiments (Figure A2c Ivins and James, 1999, Geophys. J. Int.)
 md = triangle(model(), '../Exp/RoundFrontEISMINT.exp', 200000)
 md = setmask(md, '', '')
-md = parameterize(md, '../Par/GiaIvinsBenchmarksAB.py')
+md = parameterize(md, '../Par/GiaIvinsBenchmarksCD.py')
 
 #indicate what you want to compute
-md.gia.cross_section_shape = 1  # for square-edged x - section
+md.gia.cross_section_shape = 1  # for square-edged x-section
 
 #define loading history
-md.timestepping.start_time = 2005100  # after 5 kyr of deglaciation
-md.timestepping.final_time = 2500000  # 2, 500 kyr
+md.timestepping.start_time = 0.3  # for t \approx 0 kyr : to get eleastic response!
+md.timestepping.final_time = 2500000  # 2,500 kyr
 md.geometry.thickness = np.array([np.append(md.geometry.thickness * 0.0, 0.0),
-                                  np.append(md.geometry.thickness, 1000),
-                                  np.append(md.geometry.thickness, 2000000),
-                                  np.append(md.geometry.thickness * 0.0, 2000100),
-                                  np.append(md.geometry.thickness * 0.0, md.timestepping.start_time)]).T
+                                  np.append(md.geometry.thickness / 2.0, 0.1),
+                                  np.append(md.geometry.thickness, 0.2),
+                                  np.append(md.geometry.thickness, md.timestepping.start_time)]).T
 
 #solve for GIA deflection
@@ -34,5 +33,5 @@
 
 #Fields and tolerances to track changes
-field_names = ['Ugia', 'UGiaRate']
+field_names = ['UGia', 'UGiaRate']
 field_tolerances = [1e-13, 1e-13]
 field_values = [md.results.GiaSolution.UGia, md.results.GiaSolution.UGiaRate]
Index: /issm/trunk-jpl/test/NightlyRun/test2053.m
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test2053.m	(revision 25289)
+++ /issm/trunk-jpl/test/NightlyRun/test2053.m	(revision 25290)
@@ -1,32 +1,63 @@
-%Test Name: GiaIvinsBenchmarksAB2dA3
+%Test Name: GiaIvinsBenchmarksAB2dD
 % Benchmark experiments (Figure A2a Ivins and James, 1999, Geophys. J. Int.) 
 md=triangle(model(),'../Exp/RoundFrontEISMINT.exp',200000.);
 md=setmask(md,'','');
-md=parameterize(md,'../Par/GiaIvinsBenchmarksAB.par');
+md=parameterize(md,'../Par/GiaIvinsBenchmarksCD.par');
 
 %% indicate what you want to compute 
-md.gia.cross_section_shape=1;    % for square-edged x-section 
+md.gia.cross_section_shape=2;    % for elliptical edge  
 
-%% define loading history 
-md.timestepping.start_time=2010100; % after 10 kyr of deglaciation 
+% evaluation time (termed start_time) 
+md.timestepping.start_time=0.3;     % for t \approx 0 kyr : to get eleastic response!   
 md.timestepping.final_time=2500000; % 2,500 kyr
+
+%% define loading history {{{ 
 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],...
+	[md.geometry.thickness/2.0; 0.1],...
+	[md.geometry.thickness; 0.2],...
+	[md.geometry.thickness; md.timestepping.start_time],...
 	];
+% }}} 
+
+% find out elements that have zero loads throughout the loading history.
+pos = find(sum(abs(md.geometry.thickness(1:end-1,:)),2)==0);
+md.mask.ice_levelset(pos)=1; % no-ice.
+
+md.cluster=generic('name',oshostname(),'np',3);
+md.verbose=verbose('1111111');
 
 %% solve for GIA deflection 
-md.cluster=generic('name',oshostname(),'np',3);
-md.verbose=verbose('1111111');
 md=solve(md,'Gia');
 
+%Test Name: GiaIvinsBenchmarksAB2dD1
+U_AB2dD1 = md.results.GiaSolution.UGia; 
+URate_AB2dD1 = md.results.GiaSolution.UGiaRate; 
+
+%Test Name: GiaIvinsBenchmarksAB2dD2
+%% different evaluation time. {{{ 
+md.timestepping.start_time=1000.3;  % for t \approx 1 kyr 
+md.geometry.thickness(end,end) = md.timestepping.start_time; 
+
+md=solve(md,'Gia');
+
+U_AB2dD2 = md.results.GiaSolution.UGia; 
+URate_AB2dD2 = md.results.GiaSolution.UGiaRate; 
+% }}} 
+
+%Test Name: GiaIvinsBenchmarksAB2dD3
+%% different evaluation time. {{{ 
+md.timestepping.start_time=2400000; % for t \approx \infty 
+md.geometry.thickness(end,end) = md.timestepping.start_time; 
+
+md=solve(md,'Gia');
+
+U_AB2dD3 = md.results.GiaSolution.UGia; 
+URate_AB2dD3 = md.results.GiaSolution.UGiaRate; 
+% }}} 
+
 %Fields and tolerances to track changes
-field_names     ={'UGia','UGiaRate'};
-field_tolerances={1e-13,1e-13};
-field_values={...
-	(md.results.GiaSolution.UGia),...
-	(md.results.GiaSolution.UGiaRate),...
-	};
+field_names     ={'U_AB2dD1','URate_AB2dD1','U_AB2dD2','URate_AB2dD2','U_AB2dD3','URate_AB2dD3'};
+field_tolerances={1e-13,1e-13,1e-13,1e-13,1e-13,1e-13};
+field_values={U_AB2dD1,URate_AB2dD1,U_AB2dD2,URate_AB2dD2,U_AB2dD3,URate_AB2dD3}; 
+
Index: /issm/trunk-jpl/test/NightlyRun/test2053.py
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test2053.py	(revision 25289)
+++ /issm/trunk-jpl/test/NightlyRun/test2053.py	(revision 25290)
@@ -1,3 +1,3 @@
-#Test Name: GiaIvinsBenchmarksAB2dA3
+#Test Name: GiaIvinsBenchmarksAB2dD1
 from socket import gethostname
 
@@ -11,20 +11,19 @@
 
 
-#Benchmark experiments (Figure A2a Ivins and James, 1999, Geophys. J. Int.)
+#Benchmark experiments (Figure A2c Ivins and James, 1999, Geophys. J. Int.)
 md = triangle(model(), '../Exp/RoundFrontEISMINT.exp', 200000)
 md = setmask(md, '', '')
-md = parameterize(md, '../Par/GiaIvinsBenchmarksAB.py')
+md = parameterize(md, '../Par/GiaIvinsBenchmarksCD.py')
 
 #indicate what you want to compute
-md.gia.cross_section_shape = 1  # for square-edged x - section
+md.gia.cross_section_shape = 2  # for square-edged x - section
 
 #define loading history
-md.timestepping.start_time = 2010100  # after 10 kyr of deglaciation
+md.timestepping.start_time = 0.3  # for t \approx 0 kyr : to get eleastic response!
 md.timestepping.final_time = 2500000  # 2, 500 kyr
 md.geometry.thickness = np.array([np.append(md.geometry.thickness * 0.0, 0.0),
-                                  np.append(md.geometry.thickness, 1000),
-                                  np.append(md.geometry.thickness, 2000000),
-                                  np.append(md.geometry.thickness * 0.0, 2000100),
-                                  np.append(md.geometry.thickness * 0.0, md.timestepping.start_time)]).T
+                                  np.append(md.geometry.thickness / 2.0, 0.1),
+                                  np.append(md.geometry.thickness, 0.2),
+                                  np.append(md.geometry.thickness, md.timestepping.start_time)]).T
 
 #solve for GIA deflection
@@ -34,5 +33,5 @@
 
 #Fields and tolerances to track changes
-field_names = ['Ugia', 'UgiaRate']
+field_names = ['UGia', 'UGiaRate']
 field_tolerances = [1e-13, 1e-13]
 field_values = [md.results.GiaSolution.UGia, md.results.GiaSolution.UGiaRate]
Index: /issm/trunk-jpl/test/Par/GiaIvinsBenchmarksAB.par
===================================================================
--- /issm/trunk-jpl/test/Par/GiaIvinsBenchmarksAB.par	(revision 25289)
+++ /issm/trunk-jpl/test/Par/GiaIvinsBenchmarksAB.par	(revision 25290)
@@ -7,5 +7,5 @@
 		md.geometry.thickness(i)=2000.0;
 	else
-		md.geometry.thickness(i)=1.0; % non-zero thickness
+		md.geometry.thickness(i)=0; 
 	end
 end
Index: /issm/trunk-jpl/test/Par/GiaIvinsBenchmarksCD.par
===================================================================
--- /issm/trunk-jpl/test/Par/GiaIvinsBenchmarksCD.par	(revision 25289)
+++ /issm/trunk-jpl/test/Par/GiaIvinsBenchmarksCD.par	(revision 25290)
@@ -7,5 +7,5 @@
 		md.geometry.thickness(i)=3000.0;
 	else
-		md.geometry.thickness(i)=1.0; % non-zero thickness
+		md.geometry.thickness(i)=0; 
 	end
 end