Index: /issm/trunk-jpl/src/m/boundaryconditions/getlovenumbers.py
===================================================================
--- /issm/trunk-jpl/src/m/boundaryconditions/getlovenumbers.py	(revision 25182)
+++ /issm/trunk-jpl/src/m/boundaryconditions/getlovenumbers.py	(revision 25183)
@@ -10069,9 +10069,9 @@
     elif frame == 'CF': #from Blewitt, 2003, JGR
         if type == 'loadingverticaldisplacement':
-            series[0] = -0.269
+            series[1] = -0.269
         elif type == 'loadinggravitationalpotential':
-            series[0] = 0.021
+            series[1] = 0.021
         elif type == 'loadinghorizontaldisplacement':
-            series[0] = 0.134
+            series[1] = 0.134
     else:
         raise Exception("love_numbers error message: unknown reference frame: {}".format(frame))
Index: /issm/trunk-jpl/src/m/classes/lovenumbers.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/lovenumbers.py	(revision 25182)
+++ /issm/trunk-jpl/src/m/classes/lovenumbers.py	(revision 25183)
@@ -78,4 +78,6 @@
         s += '{}\n'.format(fielddisplay(self, 'tk', 'tidal load Love number (deg 2)'))
         s += '{}\n'.format(fielddisplay(self, 'tk2secular', 'secular fluid Love number'))
+
+        return s
     #}}}
 
Index: /issm/trunk-jpl/test/NightlyRun/test2111.m
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test2111.m	(revision 25182)
+++ /issm/trunk-jpl/test/NightlyRun/test2111.m	(revision 25183)
@@ -3,15 +3,16 @@
 
 %mesh ais: {{{
-	md=model();
-	md=triangle(md,'../Exp/Ais.exp',200000); % max element size
+md=model();
+md=triangle(md,'../Exp/Ais.exp',200000); % max element size
 % }}}
 %define load: {{{
-	md.esa.deltathickness=zeros(md.mesh.numberofelements,1);
-	disc_radius=500; % km
-	index=md.mesh.elements;
-	x_element=mean(md.mesh.x(index),2)-1.0e6;
-	y_element=mean(md.mesh.y(index),2)-1.0e6;
-	rad_dist=sqrt(x_element.^2+y_element.^2)/1000;  % radial distance in km
-	md.esa.deltathickness(rad_dist<=disc_radius)=-1;   % 1 m water withdrawl
+md.esa.deltathickness=zeros(md.mesh.numberofelements,1);
+disc_radius=500; % km
+index=md.mesh.elements;
+x_element=mean(md.mesh.x(index),2)-1.0e6;
+y_element=mean(md.mesh.y(index),2)-1.0e6;
+rad_dist=sqrt(x_element.^2+y_element.^2)/1000;  % radial distance in km
+pos=find(rad_dist<=disc_radius);
+md.esa.deltathickness(pos)=-1;   % 1 m water withdrawl
 % }}}
 %read in love numbers:{{{
@@ -19,28 +20,30 @@
 % }}}
 %mask:  {{{
-	%make sure wherever there is an ice load, that the mask is set to ice: 
-	md.mask.ice_levelset=ones(md.mesh.numberofvertices,1);
-	pos=find(md.esa.deltathickness); md.mask.ice_levelset(md.mesh.elements(pos,:))=-1;
+%make sure wherever there is an ice load, that the mask is set to ice: 
+md.mask.ice_levelset=ones(md.mesh.numberofvertices,1);
+pos=find(md.esa.deltathickness);
+md.mask.ice_levelset(md.mesh.elements(pos,:))=-1;
 
-	%is ice grounded? 
-	md.mask.ocean_levelset=-ones(md.mesh.numberofvertices,1);
-	pos=find(md.mask.ice_levelset<=0); md.mask.ocean_levelset(pos)=1; 
+%is ice grounded? 
+md.mask.ocean_levelset=-ones(md.mesh.numberofvertices,1);
+pos=find(md.mask.ice_levelset<=0);
+md.mask.ocean_levelset(pos)=1;
 % }}}
 %geometry:  {{{
-	di=md.materials.rho_ice/md.materials.rho_water;
-	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;
+di=md.materials.rho_ice/md.materials.rho_water;
+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);
+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);
 % }}}
 %additional parameters, miscellaneous: {{{
-	md.miscellaneous.name='test2111';
-	md.esa.degacc=0.01;
-	md.esa.hemisphere = -1; 
+md.miscellaneous.name='test2111';
+md.esa.degacc=0.01;
+md.esa.hemisphere = -1; 
 % }}}
 
Index: /issm/trunk-jpl/test/NightlyRun/test2111.py
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test2111.py	(revision 25182)
+++ /issm/trunk-jpl/test/NightlyRun/test2111.py	(revision 25183)
@@ -17,11 +17,12 @@
 # }}}
 #define load: {{{
-md.esa.deltathickness = np.zeros((md.mesh.numberofelements, ))
+md.esa.deltathickness = np.zeros((md.mesh.numberofelements, 1))
 disc_radius = 500  # km
 index = md.mesh.elements
-x_element = np.mean(md.mesh.x[index - 1], 1) - 1.0e6
-y_element = np.mean(md.mesh.y[index - 1], 1) - 1.0e6
+x_element = md.mesh.x[index - 1].mean(axis=1) - 1.0e6
+y_element = md.mesh.y[index - 1].mean(axis=1) - 1.0e6
 rad_dist = np.sqrt(x_element**2 + y_element**2) / 1000  # radial distance in km
-md.esa.deltathickness[np.where(rad_dist <= disc_radius)] = -1  # 1 m water withdrawl
+pos = np.where(rad_dist <= disc_radius)[0]
+md.esa.deltathickness[pos] = -1  # 1 m water withdrawl
 # }}}
 #love numbers: {{{
@@ -30,24 +31,24 @@
 #mask:  {{{
 #make sure wherever there is an ice load, that the mask is set to ice:
-md.mask.ice_levelset = np.ones((md.mesh.numberofvertices, ))
-pos = np.where(md.esa.deltathickness)
+md.mask.ice_levelset = np.ones((md.mesh.numberofvertices, 1))
+pos = np.where(md.esa.deltathickness)[0]
 md.mask.ice_levelset[md.mesh.elements[pos, :]] = -1
 
 #is ice grounded?
-md.mask.ocean_levelset = -np.ones((md.mesh.numberofvertices, ))
-pos = np.where(md.mask.ice_levelset <= 0)
+md.mask.ocean_levelset = -np.ones((md.mesh.numberofvertices, 1))
+pos = np.where(md.mask.ice_levelset <= 0)[0]
 md.mask.ocean_levelset[pos] = 1
 # }}}
 #geometry:  {{{
 di = md.materials.rho_ice / md.materials.rho_water
-md.geometry.thickness = np.ones((md.mesh.numberofvertices, ))
-md.geometry.surface = (1 - di) * np.zeros((md.mesh.numberofvertices, ))
+md.geometry.thickness = np.ones((md.mesh.numberofvertices, 1))
+md.geometry.surface = (1 - di) * np.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 * np.ones((md.mesh.numberofvertices, ))
+md.initialization.temperature = 273.25 * np.ones((md.mesh.numberofvertices, 1))
 md.materials.rheology_B = paterson(md.initialization.temperature)
-md.materials.rheology_n = 3 * np.ones((md.mesh.numberofelements, ))
+md.materials.rheology_n = 3 * np.ones((md.mesh.numberofelements, 1))
 # }}}
 #additional parameters, miscellaneous: {{{
Index: /issm/trunk-jpl/test/NightlyRun/test2112.m
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test2112.m	(revision 25182)
+++ /issm/trunk-jpl/test/NightlyRun/test2112.m	(revision 25183)
@@ -3,15 +3,16 @@
 
 %mesh ais: {{{
-	md=model();
-	md=triangle(md,'../Exp/Ais.exp',100000); % max element size
+md=model();
+md=triangle(md,'../Exp/Ais.exp',100000); % max element size
 % }}}
 %define load: {{{
-	md.esa.deltathickness=zeros(md.mesh.numberofelements,1);
-	disc_radius=500; % km
-	index=md.mesh.elements;
-	x_element=mean(md.mesh.x(index),2); 
-	y_element=mean(md.mesh.y(index),2); 
-	rad_dist=sqrt(x_element.^2+y_element.^2)/1000;  % radial distance in km
-	md.esa.deltathickness(rad_dist<=disc_radius)=-1;   % 1 m water withdrawl
+md.esa.deltathickness=zeros(md.mesh.numberofelements,1);
+disc_radius=500; % km
+index=md.mesh.elements;
+x_element=mean(md.mesh.x(index),2); 
+y_element=mean(md.mesh.y(index),2); 
+rad_dist=sqrt(x_element.^2+y_element.^2)/1000;  % radial distance in km
+pos=find(rad_dist<=disc_radius);
+md.esa.deltathickness(pos)=-1;   % 1 m water withdrawl
 % }}}
 %read in love numbers:{{{
@@ -19,28 +20,28 @@
 % }}}
 %mask:  {{{
-	%make sure wherever there is an ice load, that the mask is set to ice: 
-	md.mask.ice_levelset=ones(md.mesh.numberofvertices,1);
-	pos=find(md.esa.deltathickness); md.mask.ice_levelset(md.mesh.elements(pos,:))=-1;
+%make sure wherever there is an ice load, that the mask is set to ice: 
+md.mask.ice_levelset=ones(md.mesh.numberofvertices,1);
+pos=find(md.esa.deltathickness); md.mask.ice_levelset(md.mesh.elements(pos,:))=-1;
 
-	%is ice grounded? 
-	md.mask.ocean_levelset=-ones(md.mesh.numberofvertices,1);
-	pos=find(md.mask.ice_levelset<=0); md.mask.ocean_levelset(pos)=1; 
+%is ice grounded? 
+md.mask.ocean_levelset=-ones(md.mesh.numberofvertices,1);
+pos=find(md.mask.ice_levelset<=0); md.mask.ocean_levelset(pos)=1; 
 % }}}
 %geometry:  {{{
-	di=md.materials.rho_ice/md.materials.rho_water;
-	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;
+di=md.materials.rho_ice/md.materials.rho_water;
+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);
+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);
 % }}}
 %additional parameters, miscellaneous: {{{
-	md.miscellaneous.name='test2112';
-	md.esa.degacc=0.01;
-	md.esa.hemisphere = -1; 
+md.miscellaneous.name='test2112';
+md.esa.degacc=0.01;
+md.esa.hemisphere = -1; 
 % }}}
 
Index: /issm/trunk-jpl/test/NightlyRun/test2112.py
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test2112.py	(revision 25182)
+++ /issm/trunk-jpl/test/NightlyRun/test2112.py	(revision 25183)
@@ -20,8 +20,9 @@
 disc_radius = 500  # km
 index = md.mesh.elements
-x_element = np.mean(md.mesh.x[index - 1], 1)
-y_element = np.mean(md.mesh.y[index - 1], 1)
+x_element = md.mesh.x[index - 1].mean(axis=1)
+y_element = md.mesh.y[index - 1].mean(axis=1)
 rad_dist = np.sqrt(x_element**2 + y_element**2) / 1000  # radial distance in km
-md.esa.deltathickness[np.where(rad_dist <= disc_radius)] = -1  # 1 m water withdrawl
+pos = np.where(rad_dist <= disc_radius)[0]
+md.esa.deltathickness[pos] = -1  # 1 m water withdrawl
 # }}}
 #love numbers: {{{
@@ -30,24 +31,24 @@
 #mask:  {{{
 #make sure wherever there is an ice load, that the mask is set to ice:
-md.mask.ice_levelset = np.ones((md.mesh.numberofvertices, ))
-pos = np.where(md.esa.deltathickness)
+md.mask.ice_levelset = np.ones((md.mesh.numberofvertices, 1))
+pos = np.nonzero(md.esa.deltathickness)[0]
 md.mask.ice_levelset[md.mesh.elements[pos, :]] = -1
 
 #is ice grounded?
-md.mask.ocean_levelset = -np.ones((md.mesh.numberofvertices, ))
-pos = np.where(md.mask.ice_levelset <= 0)
+md.mask.ocean_levelset = -np.ones((md.mesh.numberofvertices, 1))
+pos = np.where(md.mask.ice_levelset <= 0)[0]
 md.mask.ocean_levelset[pos] = 1
 # }}}
 #geometry:  {{{
 di = md.materials.rho_ice / md.materials.rho_water
-md.geometry.thickness = np.ones((md.mesh.numberofvertices, ))
-md.geometry.surface = (1 - di) * np.zeros((md.mesh.numberofvertices, ))
+md.geometry.thickness = np.ones((md.mesh.numberofvertices, 1))
+md.geometry.surface = (1 - di) * np.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 * np.ones((md.mesh.numberofvertices, ))
+md.initialization.temperature = 273.25 * np.ones((md.mesh.numberofvertices, 1))
 md.materials.rheology_B = paterson(md.initialization.temperature)
-md.materials.rheology_n = 3 * np.ones((md.mesh.numberofelements, ))
+md.materials.rheology_n = 3 * np.ones((md.mesh.numberofelements, 1))
 # }}}
 #additional parameters, miscellaneous: {{{
Index: /issm/trunk-jpl/test/NightlyRun/test2113.m
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test2113.m	(revision 25182)
+++ /issm/trunk-jpl/test/NightlyRun/test2113.m	(revision 25183)
@@ -4,15 +4,16 @@
 
 %mesh ais: {{{
-	md=model();
-	md=triangle(md,'../Exp/Ais.exp',200000); % max element size
+md=model();
+md=triangle(md,'../Exp/Ais.exp',200000); % max element size
 % }}}
 %define load: {{{
-	md.esa.deltathickness=zeros(md.mesh.numberofelements,1);
-	disc_radius=500; % km
-	index=md.mesh.elements;
-	x_element=mean(md.mesh.x(index),2)-1.0e6;
-	y_element=mean(md.mesh.y(index),2)+1.0e6;
-	rad_dist=sqrt(x_element.^2+y_element.^2)/1000;  % radial distance in km
-	md.esa.deltathickness(rad_dist<=disc_radius)=-1;   % 1 m water withdrawl
+md.esa.deltathickness=zeros(md.mesh.numberofelements,1);
+disc_radius=500; % km
+index=md.mesh.elements;
+x_element=mean(md.mesh.x(index),2)-1.0e6;
+y_element=mean(md.mesh.y(index),2)+1.0e6;
+rad_dist=sqrt(x_element.^2+y_element.^2)/1000;  % radial distance in km
+pos=find(rad_dist<=disc_radius);
+md.esa.deltathickness(pos)=-1;   % 1 m water withdrawl
 % }}}
 %read in love numbers:{{{
@@ -20,28 +21,28 @@
 % }}}
 %mask:  {{{
-	%make sure wherever there is an ice load, that the mask is set to ice: 
-	md.mask.ice_levelset=ones(md.mesh.numberofvertices,1);
-	pos=find(md.esa.deltathickness); md.mask.ice_levelset(md.mesh.elements(pos,:))=-1;
+%make sure wherever there is an ice load, that the mask is set to ice: 
+md.mask.ice_levelset=ones(md.mesh.numberofvertices,1);
+pos=find(md.esa.deltathickness); md.mask.ice_levelset(md.mesh.elements(pos,:))=-1;
 
-	%is ice grounded? 
-	md.mask.ocean_levelset=-ones(md.mesh.numberofvertices,1);
-	pos=find(md.mask.ice_levelset<=0); md.mask.ocean_levelset(pos)=1; 
+%is ice grounded? 
+md.mask.ocean_levelset=-ones(md.mesh.numberofvertices,1);
+pos=find(md.mask.ice_levelset<=0); md.mask.ocean_levelset(pos)=1; 
 % }}}
 %geometry:  {{{
-	di=md.materials.rho_ice/md.materials.rho_water;
-	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;
+di=md.materials.rho_ice/md.materials.rho_water;
+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);
+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);
 % }}}
 %additional parameters, miscellaneous: {{{
-	md.miscellaneous.name='test2113';
-	md.esa.degacc=0.01;
-	md.esa.hemisphere = 1; % AIS is placed in Northern Hemisphere 
+md.miscellaneous.name='test2113';
+md.esa.degacc=0.01;
+md.esa.hemisphere = 1; % AIS is placed in Northern Hemisphere 
 % }}}
 
Index: /issm/trunk-jpl/test/NightlyRun/test2113.py
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test2113.py	(revision 25182)
+++ /issm/trunk-jpl/test/NightlyRun/test2113.py	(revision 25183)
@@ -18,5 +18,5 @@
 # }}}
 #define load: {{{
-md.esa.deltathickness = np.zeros((md.mesh.numberofelements, ))
+md.esa.deltathickness = np.zeros((md.mesh.numberofelements, 1))
 disc_radius = 500  # km
 index = md.mesh.elements
@@ -24,5 +24,6 @@
 y_element = np.mean(md.mesh.y[index - 1], 1) + 1.0e6
 rad_dist = np.sqrt(x_element**2 + y_element**2) / 1000  # radial distance in km
-md.esa.deltathickness[np.where(rad_dist <= disc_radius)] = -1  # 1 m water withdrawl
+pos = np.where(rad_dist <= disc_radius)[0]
+md.esa.deltathickness[pos] = -1  # 1 m water withdrawl
 # }}}
 #love numbers: {{{
@@ -31,24 +32,24 @@
 #mask:  {{{
 #make sure wherever there is an ice load, that the mask is set to ice:
-md.mask.ice_levelset = np.ones((md.mesh.numberofvertices, ))
-pos = np.where(md.esa.deltathickness)
+md.mask.ice_levelset = np.ones((md.mesh.numberofvertices, 1))
+pos = np.nonzero(md.esa.deltathickness)[0]
 md.mask.ice_levelset[md.mesh.elements[pos, :]] = -1
 
 #is ice grounded?
-md.mask.ocean_levelset = -np.ones((md.mesh.numberofvertices, ))
-pos = np.where(md.mask.ice_levelset <= 0)
+md.mask.ocean_levelset = -np.ones((md.mesh.numberofvertices, 1))
+pos = np.where(md.mask.ice_levelset <= 0)[0]
 md.mask.ocean_levelset[pos] = 1
 # }}}
 #geometry:  {{{
 di = md.materials.rho_ice / md.materials.rho_water
-md.geometry.thickness = np.ones((md.mesh.numberofvertices, ))
-md.geometry.surface = (1 - di) * np.zeros((md.mesh.numberofvertices, ))
+md.geometry.thickness = np.ones((md.mesh.numberofvertices, 1))
+md.geometry.surface = (1 - di) * np.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 * np.ones((md.mesh.numberofvertices, ))
+md.initialization.temperature = 273.25 * np.ones((md.mesh.numberofvertices, 1))
 md.materials.rheology_B = paterson(md.initialization.temperature)
-md.materials.rheology_n = 3 * np.ones((md.mesh.numberofelements, ))
+md.materials.rheology_n = 3 * np.ones((md.mesh.numberofelements, 1))
 # }}}
 #additional parameters, miscellaneous: {{{