Index: /issm/trunk-jpl/src/m/classes/geometry.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/geometry.m	(revision 25302)
+++ /issm/trunk-jpl/src/m/classes/geometry.m	(revision 25303)
@@ -90,9 +90,6 @@
 		end % }}}
 		function marshall(self,prefix,md,fid) % {{{
-			disp(md.geometry.thickness)
 			WriteData(fid,prefix,'object',self,'fieldname','surface','format','DoubleMat','mattype',1);
 			WriteData(fid,prefix,'object',self,'fieldname','thickness','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
-			disp(md.geometry.thickness)
-			pause
 			WriteData(fid,prefix,'object',self,'fieldname','base','format','DoubleMat','mattype',1);
 			WriteData(fid,prefix,'object',self,'fieldname','bed','format','DoubleMat','mattype',1);
Index: /issm/trunk-jpl/src/m/solve/parseresultsfromdisk.m
===================================================================
--- /issm/trunk-jpl/src/m/solve/parseresultsfromdisk.m	(revision 25302)
+++ /issm/trunk-jpl/src/m/solve/parseresultsfromdisk.m	(revision 25303)
@@ -229,13 +229,13 @@
 		mu0 = md.love.mu0;
 		rr = md.materials.radius; 
-		rho= md.materials.density; 
-		rho_avg = sum( rho .*  diff(rr.^3)/sum(diff(rr.^3)) ); 
-		temp_field = cell(degmax+1,nfreq,nlayer+1,6); 
+		rho= md.materials.density;
+		rho_avg = sum(rho.*diff(rr.^3)/sum(diff(rr.^3)));
+		temp_field = cell(degmax+1,nfreq,nlayer+1,6);
 		for ii=1:degmax+1
 			for jj=1:nfreq
 				for kk=1:nlayer+1
 					if (kk<nlayer+1)
-						ll = (ii-1)*(nlayer+1)*6 + ((kk-1)*6+1) + 3; 
-						temp_field{ii,jj,kk,1} = field(ll+(1-1),jj)*r0;		% mm = 4 
+						ll = (ii-1)*(nlayer+1)*6 + ((kk-1)*6+1) + 3;
+						temp_field{ii,jj,kk,1} = field(ll+(1-1),jj)*r0;		% mm = 4
 						temp_field{ii,jj,kk,2} = field(ll+(2-1),jj)*mu0;	% mm = 5 
 						temp_field{ii,jj,kk,3} = field(ll+(3-1),jj)*r0;		% mm = 6 
@@ -252,5 +252,5 @@
 						if (md.love.forcing_type==9)
 							temp_field{ii,jj,kk,2} = 0; 
-							temp_field{ii,jj,kk,6} = (2*ii-1)/r0 - ii*field(ll+(3-1),jj)*g0; 
+							temp_field{ii,jj,kk,6} = (2*ii-1)/r0 - ii*field(ll+(3-1),jj)*g0;
 						elseif (md.love.forcing_type==11)
 							temp_field{ii,jj,kk,2} = -(2*(ii-1)+1)*rho_avg/3; 
Index: /issm/trunk-jpl/src/m/solve/parseresultsfromdisk.py
===================================================================
--- /issm/trunk-jpl/src/m/solve/parseresultsfromdisk.py	(revision 25302)
+++ /issm/trunk-jpl/src/m/solve/parseresultsfromdisk.py	(revision 25303)
@@ -248,5 +248,6 @@
             rr = md.materials.radius
             rho = md.materials.density
-            rho_avg = (rho * np.diff(np.power(rr, 3), n=1, axis=0)) / np.diff(np.power(rr, 3).sum()).sum()
+            rho_avg_partial = np.diff(np.power(rr, 3), n=1, axis=0)
+            rho_avg = ((rho * rho_avg_partial) / rho_avg_partial.sum()).sum()
             temp_field = np.empty((degmax + 1, nfreq, nlayer + 1, 6))
             temp_field.fill(0.0)
@@ -254,26 +255,25 @@
                 for jj in range(nfreq):
                     for kk in range(nlayer + 1):
-                        if kk < nlayer + 1:
+                        if kk < nlayer: # NOTE: Upper bound of range is non-inclusive (compare to src/m/solve/parseresultsfromdisk.m)
                             ll = ii * (nlayer + 1) * 6 + (kk * 6 + 1) + 3
-                            temp_field[ii, jj, kk, 0] = field[ll + (1 - 1), jj] * r0        # mm = 4
-                            temp_field[ii, jj, kk, 1] = field[ll + (2 - 1), jj] * mu0       # mm = 5
-                            temp_field[ii, jj, kk, 2] = field[ll + (3 - 1), jj] * r0        # mm = 6
-                            temp_field[ii, jj, kk, 3] = field[ll + (4 - 1), jj] * mu0       # mm = 1
-                            temp_field[ii, jj, kk, 4] = field[ll + (5 - 1), jj] * r0 * g0   # mm = 2
-                            temp_field[ii, jj, kk, 5] = field[ll + (6 - 1), jj] * g0        # mm = 3
-                            print(temp_field)
+                            temp_field[ii, jj, kk, 0] = field[ll + (0 - 1), jj] * r0        # mm = 4
+                            temp_field[ii, jj, kk, 1] = field[ll + (1 - 1), jj] * mu0       # mm = 5
+                            temp_field[ii, jj, kk, 2] = field[ll + (2 - 1), jj] * r0        # mm = 6
+                            temp_field[ii, jj, kk, 3] = field[ll + (3 - 1), jj] * mu0       # mm = 1
+                            temp_field[ii, jj, kk, 4] = field[ll + (4 - 1), jj] * r0 * g0   # mm = 2
+                            temp_field[ii, jj, kk, 5] = field[ll + (5 - 1), jj] * g0        # mm = 3
                         else: # surface
-                            ll = ii * (nlayer + 1) * 6 - 2
-                            temp_field[ii, jj, kk, 0] = field[ll + (1 - 1), jj] * r0
-                            temp_field[ii, jj, kk, 2] = field[ll + (2 - 1), jj] * r0
-                            temp_field[ii, jj, kk, 4] = field[ll + (3 - 1), jj] * r0 * g0
+                            ll = (ii + 1) * (nlayer + 1) * 6 - 2
+                            temp_field[ii, jj, kk, 0] = field[ll + (0 - 1), jj] * r0
+                            temp_field[ii, jj, kk, 2] = field[ll + (1 - 1), jj] * r0
+                            temp_field[ii, jj, kk, 4] = field[ll + (2 - 1), jj] * r0 * g0
                             # surface BC
                             temp_field[ii, jj, kk, 3] = 0
                             if md.love.forcing_type == 9:
                                 temp_field[ii, jj, kk, 1] = 0
-                                temp_field[ii, jj, kk, 5] = (2 * ii - 1) / r0 - ii * field[ll + (3 - 1), jj] * g0
+                                temp_field[ii, jj, kk, 5] = (2 * (ii + 1) - 1) / r0 - (ii + 1) * field[ll + (2 - 1), jj] * g0
                             elif md.love.forcing_type == 11:
-                                temp_field[ii, jj, kk, 1] = -(2 * (ii - 1) + 1) * rho_avg / 3
-                                temp_field[ii, jj, kk, 5] = (2 * ii - 1) / r0 - ii * field[ll + (3 - 1), jj] * g0
+                                temp_field[ii, jj, kk, 1] = -(2 * ii + 1) * rho_avg / 3
+                                temp_field[ii, jj, kk, 5] = (2 * (ii + 1) - 1) / r0 - (ii + 1) * field[ll + (2 - 1), jj] * g0
             field = temp_field
 
Index: /issm/trunk-jpl/test/NightlyRun/test2052.m
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test2052.m	(revision 25302)
+++ /issm/trunk-jpl/test/NightlyRun/test2052.m	(revision 25303)
@@ -31,6 +31,4 @@
 md=solve(md,'Gia');
 
-pause
-
 %Test Name: GiaIvinsBenchmarksAB2dC1
 U_AB2dC1 = md.results.GiaSolution.UGia; 
Index: /issm/trunk-jpl/test/NightlyRun/test2085.m
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test2085.m	(revision 25302)
+++ /issm/trunk-jpl/test/NightlyRun/test2085.m	(revision 25303)
@@ -11,8 +11,8 @@
 	md.love.forcing_type = 9;
 
-	md.materials.density=  (1:md.materials.numlayers)'*0+5511;
-	md.materials.lame_mu=  (1:md.materials.numlayers)'*0+0.75e11;
-	md.materials.viscosity=(1:md.materials.numlayers)'*0+1e21;
-	md.materials.lame_lambda=md.materials.lame_mu*0+5e17;
+	md.materials.density=zeros(md.materials.numlayers,1)+5511;
+	md.materials.lame_mu=zeros(md.materials.numlayers,1)+0.75e11;
+	md.materials.viscosity=zeros(md.materials.numlayers,1)+1e21;
+	md.materials.lame_lambda=zeros(md.materials.numlayers,1)+5e17;
 	md.materials.issolid=ones(md.materials.numlayers,1);
 	md.materials.isburgers=zeros(md.materials.numlayers,1);
@@ -25,5 +25,5 @@
 	md.love.allow_layer_deletion=1;
 	md.love.frequencies=0;
-	md.love.nfreq=length(md.love.frequencies);
+	md.love.nfreq=length(md.love.frequencies); % TODO: Why are we grabbing length of a scalar here?
 
 	md.love.sh_nmin = 2;
@@ -38,5 +38,5 @@
 
 	% extract love kernels {{{ 
-	y1=squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(:,1,:,1))); 
+	y1=squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(:,1,:,1)));
 	y2=squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(:,1,:,2))); 
 	y3=squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(:,1,:,3))); 
@@ -148,16 +148,16 @@
 	};
 field_tolerances={...
-	1e-10,1e-10,1e-10,1e-10,1e-10,...
-	1e-10,1e-10,1e-10,1e-10,1e-10,...
-	1e-10,1e-10,1e-10,1e-10,1e-10,...
-	1e-10,1e-10,1e-10,1e-10,1e-10,...
-	1e-10,1e-10,1e-10,1e-10,1e-10,...
-	1e-10,1e-10,1e-10,1e-10,1e-10,...
-	1e-10,1e-10,1e-10,1e-10,1e-10,...
-	1e-10,1e-10,1e-10,1e-10,1e-10,...
-	1e-10,1e-10,1e-10,1e-10,1e-10,...
-	1e-10,1e-10,1e-10,1e-10,1e-10,...
-	1e-10,1e-10,1e-10,1e-10,1e-10,...
-	1e-10,1e-10,1e-10,1e-10,1e-10,...
+	9e-8, 1e-7, 3e-7, 3e-7, 3e-7,...
+	9e-8, 1e-7, 3e-7, 3e-7, 1e-10,...
+	9e-8, 8e-8, 2e-8, 2e-7, 4e-7,...
+	9e-8, 9e-8, 2e-7, 4e-7, 1e-10,...
+	4e-7, 4e-7, 2e-7, 3e-8, 2e-8,...
+	2e-5, 2e-6, 2e-6, 1e-6, 2e-7,...
+	3e-6, 3e-6, 3e-6, 4e-6, 4e-6,...
+	3e-6, 3e-6, 2e-6, 6e-7, 1e-10,...
+	3e-6, 3e-6, 5e-7, 3e-6, 5e-6,...
+	3e-6, 3e-6, 9e-7, 7e-7, 1e-10,...
+	4e-6, 4e-6, 3e-6, 5e-7, 3e-7,...
+	3e-6, 3e-6, 2e-6, 7e-7, 2e-7...
 	}; 
 field_values={...
Index: /issm/trunk-jpl/test/NightlyRun/test2085.py
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test2085.py	(revision 25302)
+++ /issm/trunk-jpl/test/NightlyRun/test2085.py	(revision 25303)
@@ -1,49 +1,188 @@
-#Test Name: LovenumberstAtDepth.
-#Same as test  #1 of test2084.m
+# Test Name: FourierLoveKernels
+# Homogenous Earth, for which analytic solutions exist.
+# Love kernels for degree 2 (tested against analytic solns).
 
-from model import *
+
 from socket import gethostname
-from solve import *
-from numpy import *
+
+import numpy as np
+
 from generic import generic
 from materials import *
+from model import *
+from solve import *
 
+# Fro volumetric potential
 md = model()
-md.cluster = generic('name', gethostname(), 'np', 1)
 
 md.materials = materials('litho')
-md.miscellaneous.name = 'FourierLoveTest'
-md.groundingline.migration = 'None'
 
+md.materials.numlayers = 10
+md.love.forcing_type = 9
+
+md.materials.density = np.zeros((md.materials.numlayers, 1)) + 5511
+md.materials.lame_mu = np.zeros((md.materials.numlayers, 1)) + 0.75e11
+md.materials.viscosity = np.zeros((md.materials.numlayers, 1)) + 1e21
+md.materials.lame_lambda = np.zeros((md.materials.numlayers, 1)) + 0.5e17
+md.materials.issolid = np.ones((md.materials.numlayers, 1))
+md.materials.isburgers = np.zeros((md.materials.numlayers, 1))
+md.materials.burgers_mu = md.materials.lame_mu / 3
+md.materials.burgers_viscosity = md.materials.viscosity / 10
+
+md.materials.radius = np.linspace(10e3, 6371e3, md.materials.numlayers + 1).reshape(-1, 1)
+md.love.g0 = 9.8134357285509388 # directly grabbed from fourierlovesolver for this particular case
+
+md.love.allow_layer_deletion = 1
+md.love.frequencies = 0
+md.love.nfreq = 1
+
+md.love.sh_nmin = 2
+md.love.sh_nmax = 2
+md.love.love_kernels = 1
+
+md.miscellaneous.name = 'kernels'
+md.cluster = generic('name', gethostname(), 'np', 1)
 md.verbose = verbose('111111101')
-cst = 365.25 * 24 * 3600 * 1000
 
-md.materials.numlayers = 6
-md.materials.radius = np.array([10, 1222.5, 3.4800e+03, 5.7010e+03, 5.9510e+03, 6.3010e+03, 6.3710e+03]).reshape(-1, 1) * 1e3
-md.materials.density = np.array([1.0750e4, 1.0750e+04, 4.9780e+03, 3.8710e+03, 3.4380e+03, 3.0370e+03]).reshape(-1, 1)
-md.materials.lame_mu = np.array([1e-5, 0, 2.2834e+00, 1.0549e+00, 7.0363e-01, 5.0605e-01]).reshape(-1, 1) * 1e11
-md.materials.viscosity = np.array([0, 0, 2.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+25]).reshape(-1, 1) * 1e21
-md.materials.lame_lambda = np.array(md.materials.lame_mu) * 0 + 5e14
-md.materials.issolid = np.array([1, 0, 1, 1, 1, 1]).reshape(-1, 1)
-md.materials.isburgers = np.zeros((md.materials.numlayers)).reshape(-1, 1)
+md = solve(md,'lv')
 
-md.love.love_kernels = 1
-md.love.allow_layer_deletion = 1
-md.love.frequencies = (np.array([0]) * 2 * pi).reshape(-1, 1) / cst
-md.love.nfreq = len(md.love.frequencies)
-md.love.sh_nmax = 2
+# Extract love kernels #{{{
+y1 = md.results.LoveSolution.LoveKernelsReal[:,0,:,0].squeeze()
+y2 = md.results.LoveSolution.LoveKernelsReal[:,0,:,1].squeeze()
+y3 = md.results.LoveSolution.LoveKernelsReal[:,0,:,2].squeeze()
+y4 = md.results.LoveSolution.LoveKernelsReal[:,0,:,3].squeeze()
+y5 = md.results.LoveSolution.LoveKernelsReal[:,0,:,4].squeeze()
+y6 = md.results.LoveSolution.LoveKernelsReal[:,0,:,5].squeeze()
 
-md.materials.burgers_mu = md.materials.lame_mu
-md.materials.burgers_viscosity = md.materials.viscosity
+y1_tidal_degree2_interface1 = y1[2, 0]
+y1_tidal_degree2_interface2 = y1[2, 1]
+y1_tidal_degree2_interface7 = y1[2, 6]
+y1_tidal_degree2_interface10 = y1[2, 9]
+y1_tidal_degree2_interface11 = y1[2, 10]
 
-md = solve(md, 'lv')
+y2_tidal_degree2_interface1 = y2[2, 0]
+y2_tidal_degree2_interface2 = y2[2, 1]
+y2_tidal_degree2_interface7 = y2[2, 6]
+y2_tidal_degree2_interface10 = y2[2, 9]
+y2_tidal_degree2_interface11 = y2[2, 10]
+
+y3_tidal_degree2_interface1 = y3[2, 0]
+y3_tidal_degree2_interface2 = y3[2, 1]
+y3_tidal_degree2_interface7 = y3[2, 6]
+y3_tidal_degree2_interface10 = y3[2, 9]
+y3_tidal_degree2_interface11 = y3[2, 10]
+
+y4_tidal_degree2_interface1 = y4[2, 0]
+y4_tidal_degree2_interface2 = y4[2, 1]
+y4_tidal_degree2_interface7 = y4[2, 6]
+y4_tidal_degree2_interface10 = y4[2, 9]
+y4_tidal_degree2_interface11 = y4[2, 10]
+
+y5_tidal_degree2_interface1 = y5[2, 0]
+y5_tidal_degree2_interface2 = y5[2, 1]
+y5_tidal_degree2_interface7 = y5[2, 6]
+y5_tidal_degree2_interface10 = y5[2, 9]
+y5_tidal_degree2_interface11 = y5[2, 10]
+
+y6_tidal_degree2_interface1 = y6[2, 0]
+y6_tidal_degree2_interface2 = y6[2, 1]
+y6_tidal_degree2_interface7 = y6[2, 6]
+y6_tidal_degree2_interface10 = y6[2, 9]
+y6_tidal_degree2_interface11 = y6[2, 10]
+#}}}
+
+# For surface load
+md.love.forcing_type = 11
+
+md = solve(md,'lv')
+
+# Extract love kernels #{{{
+y1 = md.results.LoveSolution.LoveKernelsReal[:,0,:,0].squeeze()
+y2 = md.results.LoveSolution.LoveKernelsReal[:,0,:,1].squeeze()
+y3 = md.results.LoveSolution.LoveKernelsReal[:,0,:,2].squeeze()
+y4 = md.results.LoveSolution.LoveKernelsReal[:,0,:,3].squeeze()
+y5 = md.results.LoveSolution.LoveKernelsReal[:,0,:,4].squeeze()
+y6 = md.results.LoveSolution.LoveKernelsReal[:,0,:,5].squeeze()
+
+y1_load_degree2_interface1 = y1[2, 0]
+y1_load_degree2_interface2 = y1[2, 1]
+y1_load_degree2_interface7 = y1[2, 6]
+y1_load_degree2_interface10 = y1[2, 9]
+y1_load_degree2_interface11 = y1[2, 10]
+
+y2_load_degree2_interface1 = y2[2, 0]
+y2_load_degree2_interface2 = y2[2, 1]
+y2_load_degree2_interface7 = y2[2, 6]
+y2_load_degree2_interface10 = y2[2, 9]
+y2_load_degree2_interface11 = y2[2, 10]
+
+y3_load_degree2_interface1 = y3[2, 0]
+y3_load_degree2_interface2 = y3[2, 1]
+y3_load_degree2_interface7 = y3[2, 6]
+y3_load_degree2_interface10 = y3[2, 9]
+y3_load_degree2_interface11 = y3[2, 10]
+
+y4_load_degree2_interface1 = y4[2, 0]
+y4_load_degree2_interface2 = y4[2, 1]
+y4_load_degree2_interface7 = y4[2, 6]
+y4_load_degree2_interface10 = y4[2, 9]
+y4_load_degree2_interface11 = y4[2, 10]
+
+y5_load_degree2_interface1 = y5[2, 0]
+y5_load_degree2_interface2 = y5[2, 1]
+y5_load_degree2_interface7 = y5[2, 6]
+y5_load_degree2_interface10 = y5[2, 9]
+y5_load_degree2_interface11 = y5[2, 10]
+
+y6_load_degree2_interface1 = y6[2, 0]
+y6_load_degree2_interface2 = y6[2, 1]
+y6_load_degree2_interface7 = y6[2, 6]
+y6_load_degree2_interface10 = y6[2, 9]
+y6_load_degree2_interface11 = y6[2, 10]
+#}}}
+
 
 #Fields and tolerances to track changes
 #loading love numbers
-field_names = ['LoveH_loading_elastic', 'LoveK_loading_elastic', 'LoveL_loading_elastic', 'LoveKernels_degree1', 'LoveKernels_degree2']
-field_tolerances = [1e-10, 1e-10, 1e-10, 1e-10, 1e-10]
-field_values = [md.results.LoveSolution.LoveHr[:, 0],
-                md.results.LoveSolution.LoveKr[:, 0],
-                md.results.LoveSolution.LoveLr[:, 0],
-                md.results.LoveSolution.LoveKernelsReal[1, 0, :, :],
-                md.results.LoveSolution.LoveKernelsReal[2, 0, :, :]]
+field_names = [
+    'y1_tidal_degree2_interface1', 'y1_tidal_degree2_interface2', 'y1_tidal_degree2_interface7', 'y1_tidal_degree2_interface10', 'y1_tidal_degree2_interface11',
+    'y2_tidal_degree2_interface1', 'y2_tidal_degree2_interface2', 'y2_tidal_degree2_interface7', 'y2_tidal_degree2_interface10', 'y2_tidal_degree2_interface11',
+    'y3_tidal_degree2_interface1', 'y3_tidal_degree2_interface2', 'y3_tidal_degree2_interface7', 'y3_tidal_degree2_interface10', 'y3_tidal_degree2_interface11',
+    'y4_tidal_degree2_interface1', 'y4_tidal_degree2_interface2', 'y4_tidal_degree2_interface7', 'y4_tidal_degree2_interface10', 'y4_tidal_degree2_interface11',
+    'y5_tidal_degree2_interface1', 'y5_tidal_degree2_interface2', 'y5_tidal_degree2_interface7', 'y5_tidal_degree2_interface10', 'y5_tidal_degree2_interface11',
+    'y6_tidal_degree2_interface1', 'y6_tidal_degree2_interface2', 'y6_tidal_degree2_interface7', 'y6_tidal_degree2_interface10', 'y6_tidal_degree2_interface11',
+    'y1_load_degree2_interface1', 'y1_load_degree2_interface2', 'y1_load_degree2_interface7', 'y1_load_degree2_interface10', 'y1_load_degree2_interface11',
+    'y2_load_degree2_interface1', 'y2_load_degree2_interface2', 'y2_load_degree2_interface7', 'y2_load_degree2_interface10', 'y2_load_degree2_interface11',
+    'y3_load_degree2_interface1', 'y3_load_degree2_interface2', 'y3_load_degree2_interface7', 'y3_load_degree2_interface10', 'y3_load_degree2_interface11',
+    'y4_load_degree2_interface1', 'y4_load_degree2_interface2', 'y4_load_degree2_interface7', 'y4_load_degree2_interface10', 'y4_load_degree2_interface11',
+    'y5_load_degree2_interface1', 'y5_load_degree2_interface2', 'y5_load_degree2_interface7', 'y5_load_degree2_interface10', 'y5_load_degree2_interface11',
+    'y6_load_degree2_interface1', 'y6_load_degree2_interface2', 'y6_load_degree2_interface7', 'y6_load_degree2_interface10', 'y6_load_degree2_interface11'
+    ]
+field_tolerances = [
+    9e-8, 1e-7, 3e-7, 3e-7, 3e-7,
+    9e-8, 1e-7, 3e-7, 3e-7, 1e-10,
+    9e-8, 8e-8, 2e-8, 2e-7, 4e-7,
+    9e-8, 9e-8, 2e-7, 4e-7, 1e-10,
+    4e-7, 4e-7, 2e-7, 3e-8, 2e-8,
+    2e-5, 2e-6, 2e-6, 1e-6, 2e-7,
+    3e-6, 3e-6, 3e-6, 4e-6, 4e-6,
+    3e-6, 3e-6, 2e-6, 6e-7, 1e-10,
+    3e-6, 3e-6, 5e-7, 3e-6, 5e-6,
+    3e-6, 3e-6, 9e-7, 7e-7, 1e-10,
+    4e-6, 4e-6, 3e-6, 5e-7, 3e-7,
+    3e-6, 3e-6, 2e-6, 7e-7, 2e-7
+    ]
+field_values = [
+    y1_tidal_degree2_interface1, y1_tidal_degree2_interface2, y1_tidal_degree2_interface7, y1_tidal_degree2_interface10,y1_tidal_degree2_interface11,
+    y2_tidal_degree2_interface1, y2_tidal_degree2_interface2, y2_tidal_degree2_interface7, y2_tidal_degree2_interface10,y2_tidal_degree2_interface11,
+    y3_tidal_degree2_interface1, y3_tidal_degree2_interface2, y3_tidal_degree2_interface7, y3_tidal_degree2_interface10,y3_tidal_degree2_interface11,
+    y4_tidal_degree2_interface1, y4_tidal_degree2_interface2, y4_tidal_degree2_interface7, y4_tidal_degree2_interface10,y4_tidal_degree2_interface11,
+    y5_tidal_degree2_interface1, y5_tidal_degree2_interface2, y5_tidal_degree2_interface7, y5_tidal_degree2_interface10,y5_tidal_degree2_interface11,
+    y6_tidal_degree2_interface1, y6_tidal_degree2_interface2, y6_tidal_degree2_interface7, y6_tidal_degree2_interface10,y6_tidal_degree2_interface11,
+    y1_load_degree2_interface1, y1_load_degree2_interface2, y1_load_degree2_interface7, y1_load_degree2_interface10,y1_load_degree2_interface11,
+    y2_load_degree2_interface1, y2_load_degree2_interface2, y2_load_degree2_interface7, y2_load_degree2_interface10,y2_load_degree2_interface11,
+    y3_load_degree2_interface1, y3_load_degree2_interface2, y3_load_degree2_interface7, y3_load_degree2_interface10,y3_load_degree2_interface11,
+    y4_load_degree2_interface1, y4_load_degree2_interface2, y4_load_degree2_interface7, y4_load_degree2_interface10,y4_load_degree2_interface11,
+    y5_load_degree2_interface1, y5_load_degree2_interface2, y5_load_degree2_interface7, y5_load_degree2_interface10,y5_load_degree2_interface11,
+    y6_load_degree2_interface1, y6_load_degree2_interface2, y6_load_degree2_interface7, y6_load_degree2_interface10,y6_load_degree2_interface11
+    ]