diff --git a/aux/plotting/plot_inv_results.py b/aux/plotting/plot_inv_results.py
index 75df85d..52151bd 100644
--- a/aux/plotting/plot_inv_results.py
+++ b/aux/plotting/plot_inv_results.py
@@ -39,7 +39,7 @@
 # Parameters:
 
 # Simulation Directory
-run_name = 'ismipc_rc_1e6'
+run_name = 'ismipc_1000_corr'
 dd = Path(os.environ['FENICS_ICE_BASE_DIR']) / 'example_cases' / run_name
 
 # Output Directory
@@ -167,5 +167,5 @@
 cbar.ax.set_xlabel(r'$U-U_{obs}$ (m $yr^{-1}$)')
 
 plt.tight_layout(2.0)
-plt.savefig(os.path.join(outdir, 'inv_results.pdf'))
+#plt.savefig(os.path.join(outdir, 'inv_results.pdf'))
 plt.show()
diff --git a/aux/plotting/plot_inv_results_hist.py b/aux/plotting/plot_inv_results_hist.py
index 3bc381b..b6a7351 100644
--- a/aux/plotting/plot_inv_results_hist.py
+++ b/aux/plotting/plot_inv_results_hist.py
@@ -22,6 +22,7 @@
 import matplotlib.ticker as ticker
 import seaborn as sns
 import os
+import math
 from fenics import *
 from tlm_adjoint_fenics import *
 from fenics_ice import model, config, inout
@@ -41,7 +42,7 @@
 # Parameters:
 
 # Simulation Directory
-run_name = 'ismipc_rc_1e6'
+run_name = 'ismipc_4000_corr'
 if (len(sys.argv)>1):
     run_name = sys.argv[1]
 
@@ -105,12 +106,14 @@
 uv_obs_file = str(next(results_dir.glob("*uv_obs.xml")))
 alpha_sigma_file = str(next(results_dir.glob("*alpha_sigma.xml")))
 dQ_file = str(next(results_dir.glob("*dQ_ts.xml")))
+thick_file = str(next(results_dir.glob("*H_fwd.xml")))
 
 U = Function(V, U_file)
 alpha = Function(Qp, alpha_file)
 uv_obs = Function(M, uv_obs_file)
 alpha_sigma = Function(Qp, alpha_sigma_file)
 dQ = Function(Qp,dQ_file)
+H_fwd = Function(Q,thick_file)
 # B2 = Function(M, os.path.join(dd,'B2.xml'))
 
 u, v = U.split()
@@ -167,27 +170,31 @@
 v    = B2.compute_vertex_values(mesh)
 minv = np.min(v)
 maxv = np.max(v)
-levels = np.linspace(minv,maxv,numlev)
-ticks = np.linspace(minv,maxv,3)
+levels = np.linspace(0,2000,numlev)
+ticks = np.linspace(0,2000,3)
 ax.tick_params(**tick_options)
 ax.text(0.05, 0.95, 'a', transform=ax.transAxes,
     fontsize=13, fontweight='bold', va='top')
-c = ax.tricontourf(x, y, t, v, levels = levels, cmap=plt.get_cmap(cmap))
+c = ax.tricontourf(x, y, t, v, levels = levels, cmap=plt.get_cmap(cmap),extend='both')
 cbar = plt.colorbar(c, ticks=ticks, pad=0.05, orientation="horizontal",shrink=.75)
-cbar.ax.set_xlabel(r'${\beta^2}$ (Pa $m^{-1}$ yr)')
+cbar.ax.set_xlabel(r'${C^2}$ (Pa $m^{-1}$ yr)')
 
 ax  = fig.add_subplot(232)
 ax.set_aspect('equal')
 v    = alpha_sigma.compute_vertex_values(mesh)
 minv = np.min(v)
 maxv = np.max(v)
+minv = 0
+maxv = math.ceil(np.max(v)/10.)*10.
+#levels = np.linspace(0,50,numlev)
+#ticks = np.linspace(0,50,3)
 levels = np.linspace(minv,maxv,numlev)
 ticks = np.linspace(minv,maxv,3)
 ax.tick_params(**tick_options)
 ax.text(0.05, 0.95, 'b', transform=ax.transAxes,
     fontsize=13, fontweight='bold', va='top')
 c = ax.tricontourf(x, y, v, levels = levels, cmap=plt.get_cmap(cmap))
-cbar = plt.colorbar(c, ticks=ticks, pad=0.05, orientation="horizontal", format=ticker.FormatStrFormatter('%1.1e'),shrink=.75)
+cbar = plt.colorbar(c, ticks=ticks, pad=0.15, orientation="horizontal", format=ticker.FormatStrFormatter('%1.1e'),shrink=.75)
 cbar.ax.xaxis.set_major_locator(ticker.LinearLocator(3))
 
 cbar.ax.set_xlabel(r'$\sigma$ ($Pa^{0.5}$ $m^{-0.5}$ $yr^{0.5}$)')
@@ -207,6 +214,8 @@
 
 ax  = fig.add_subplot(234)
 ax.hist(tot_err.vector()[:])
+ax.text(0.05, 0.95, 'd', transform=ax.transAxes,
+    fontsize=13, fontweight='bold', va='top')
 plt.xlabel('velocity error (m yr$^{-1}$)')
 plt.ylabel('count')
 #v   = uv_obs.compute_vertex_values(mesh)
@@ -236,6 +245,7 @@
 ax.set_aspect('equal')
 v   = dQ.compute_vertex_values(mesh)
 max_diff = np.rint(np.max(np.abs(v)))
+max_diff=5.0e8
 levels = np.linspace(-max_diff,max_diff,numlev)
 ticks = np.linspace(-max_diff,max_diff,3)
 ax.tick_params(**tick_options)
@@ -243,7 +253,22 @@
     fontsize=13, fontweight='bold', va='top')
 c = ax.tricontourf(x, y, t, v, levels = levels, cmap=plt.get_cmap(cmap_div))
 cbar = plt.colorbar(c, ticks=ticks, pad=0.05, orientation="horizontal",shrink=.75)
-cbar.ax.set_xlabel(r'$\partial Q/\partial \beta$')
+cbar.ax.set_xlabel(r'$\partial Q/\partial C$ ($m^4 Pa^{-1} m a^{-1}$)',size=10)
+
+ax  = fig.add_subplot(236)
+ax.set_aspect('equal')
+v   = H_fwd.compute_vertex_values(mesh)
+max_diff = np.rint(np.max(np.abs(1000.-v)))
+max_diff=10.
+levels = np.linspace(1000.-max_diff,1000.+max_diff,numlev)
+ticks = np.linspace(1000.-max_diff,1000.+max_diff,3)
+ax.tick_params(**tick_options)
+ax.text(0.05, 0.95, 'f', transform=ax.transAxes,
+    fontsize=13, fontweight='bold', va='top')
+c = ax.tricontourf(x, y, t, v, levels = levels, cmap=plt.get_cmap(cmap_div),extend='both')
+cbar = plt.colorbar(c, ticks=ticks, pad=0.05, orientation="horizontal",shrink=.75)
+cbar.ax.set_xlabel(r'$H(T=30)$ $(m)$')
 
 plt.tight_layout(2.0)
 plt.savefig(os.path.join(outdir, run_name + '_inv_results.png'))
+
diff --git a/aux/plotting/plotting_new/plot_eigenvalue_decay_gen_new.py b/aux/plotting/plotting_new/plot_eigenvalue_decay_gen_new.py
index 70c0d22..43d52d9 100755
--- a/aux/plotting/plotting_new/plot_eigenvalue_decay_gen_new.py
+++ b/aux/plotting/plotting_new/plot_eigenvalue_decay_gen_new.py
@@ -54,6 +54,8 @@
     lpos = np.argwhere(lam > 0)
     lneg = np.argwhere(lam < 0)
     lind = np.arange(0,len(lam))
+#    plt.semilogy(lind[lpos], lam[lpos], '.', alpha = 0.5, mew=0, label =labels[i])
+#    plt.semilogy(lind[lneg], np.abs(lam[lneg]), '.k', alpha = 0.12, mew=0,)
     plt.semilogy(lind[lpos], 1./(1+lam[lpos]), '.', alpha = 0.5, mew=0, label =labels[i])
     plt.semilogy(lind[lneg], 1./(1+np.abs(lam[lneg])), '.k', alpha = 0.12, mew=0,)
 
@@ -62,3 +64,4 @@
 plt.xlabel('Eigenvalue index')
 plt.ylabel('Uncertainty reduction')
 plt.savefig(os.path.join(outdir,figname), bbox_inches="tight")
+#plt.show()
diff --git a/aux/plotting/plotting_new/plot_paths_sample.py b/aux/plotting/plotting_new/plot_paths_sample.py
index e2186cc..5173cb7 100755
--- a/aux/plotting/plotting_new/plot_paths_sample.py
+++ b/aux/plotting/plotting_new/plot_paths_sample.py
@@ -57,7 +57,8 @@
 
     sigma_interp = np.interp(dQ_t, sigma_t, sigma_vals)
     sigma_prior_interp = np.interp(dQ_t, sigma_t, sigma_prior_vals)
-    
+   
+
     QoIs_MC = np.load(result_dir/"sampling_results.npy")
     QoIs_MC = QoIs_MC[:,QoIs_MC[-1,:]>0]
     print( str(np.shape(QoIs_MC)) + ' values OK')
@@ -74,19 +75,24 @@
     upperlim = max(upperlim,np.max(s))
     upperlim = max(upperlim,np.max(sm))
 
+    if(i==0):
+     p1 = axarr[0].plot(x, y, 'b',label='$\gamma$=50')
+     p2 = axarr[0].fill_between(x, y-s, y+s, facecolor='b', alpha = .6)
+    else:
+     p3 = axarr[0].plot(x, y, 'g',label='$\gamma$=1')
+     p4 = axarr[0].fill_between(x, y-s, y+s, facecolor='g', alpha = .6)
 
-    axarr[0].plot(x, y, 'k' if i == 0 else 'k:')
 #    axarr[0].fill_between(x, y-sp, y+sp, facecolor='b')
-    axarr[0].fill_between(x, y-s, y+s, facecolor='b' if i == 0 else 'g', alpha = .6 if i==0 else .6)
 #    axarr[0].fill_between(x, y-sm, y+sm, facecolor='r' if i == 0 else 'y', alpha = .5 if i==0 else .6)
 
-    axarr[1].semilogy(x, sp, 'b' if i == 0 else 'g')
-    axarr[1].semilogy(x, s, 'b:' if i == 0 else 'g:')
+    axarr[1].semilogy(x, sp, 'b:' if i == 0 else 'g:')
+    axarr[1].semilogy(x, s, 'b' if i == 0 else 'g')
     if (i==0):
      clr = 'blue'
     else:
      clr = 'green'
     axarr[1].semilogy(sigma_t, sigma_mc_vals,color=clr,marker='+',linestyle='none',markersize=6)
+axarr[0].legend([(p1[0],p2),(p3[0],p4)],['$\gamma$=50','$\gamma$=1'])
 
 axarr[0].set_xlabel('Time (yrs)')
 axarr[1].set_xlabel('Time (yrs)')
diff --git a/example_cases/ismipc_1000/ismipc_1000.sh b/example_cases/ismipc_1000/ismipc_1000.sh
index 5fe981d..1cde7d3 100755
--- a/example_cases/ismipc_1000/ismipc_1000.sh
+++ b/example_cases/ismipc_1000/ismipc_1000.sh
@@ -14,8 +14,8 @@ cp output_momsolve/ismipc_U_obs.h5 input/
 #Run each phase of the model in turn
 RUN_DIR=$FENICS_ICE_BASE_DIR/runs/
 python $RUN_DIR/run_inv.py ismipc_1000.toml
-#python $RUN_DIR/run_forward.py ismipc_1000.toml
-#python $RUN_DIR/run_eigendec.py ismipc_1000.toml
-#python $RUN_DIR/run_errorprop.py ismipc_1000.toml
-#python $RUN_DIR/run_invsigma.py ismipc_1000.toml
-#python $RUN_DIR/run_sample_post.py ismipc_1000.toml
+python $RUN_DIR/run_forward.py ismipc_1000.toml
+python $RUN_DIR/run_eigendec.py ismipc_1000.toml
+python $RUN_DIR/run_errorprop.py ismipc_1000.toml
+python $RUN_DIR/run_invsigma.py ismipc_1000.toml
+python $RUN_DIR/run_sample_post.py ismipc_1000.toml
diff --git a/example_cases/ismipc_1000_corr/ismipc_1000_corr.sh b/example_cases/ismipc_1000_corr/ismipc_1000_corr.sh
index 996182e..f706974 100755
--- a/example_cases/ismipc_1000_corr/ismipc_1000_corr.sh
+++ b/example_cases/ismipc_1000_corr/ismipc_1000_corr.sh
@@ -6,9 +6,9 @@ python $FENICS_ICE_BASE_DIR/aux/gen_rect_mesh.py -o ./input/momsolve_mesh.xml -x
 python $FENICS_ICE_BASE_DIR/aux/gen_rect_mesh.py -o ./input/ismip_mesh.xml -xmax 40000 -ymax 40000 -nx 30 -ny 30
 python $FENICS_ICE_BASE_DIR/aux/gen_ismipC_domain.py -o ./input/ismipc_input.h5 -L 40000 -nx 100 -ny 100 --reflect
 python $FENICS_ICE_BASE_DIR/runs/run_momsolve.py momsolve.toml
-python $FENICS_ICE_BASE_DIR/aux/Uobs_from_momsolve.py -i "U.h5" -o "ismipc_U_obs_nocorr.h5" -L 40000 -d ./output_momsolve --ls=1000. 
-cp output_momsolve/ismipc_U_obs_nocorr.h5 input/
-python $FENICS_ICE_BASE_DIR/aux/Uobs_from_momsolve_cov.py -i "U.h5" -o "ismipc_U_obs.h5" -L 40000 -d ./output_momsolve --ls=1000. -t "ismipc_1000.toml"
+#python $FENICS_ICE_BASE_DIR/aux/Uobs_from_momsolve.py -i "U.h5" -o "ismipc_U_obs_nocorr.h5" -L 40000 -d ./output_momsolve --ls=1000. 
+#cp output_momsolve/ismipc_U_obs_nocorr.h5 input/
+python $FENICS_ICE_BASE_DIR/aux/Uobs_from_momsolve.py -i "U.h5" -o "ismipc_U_obs.h5" -L 40000 -d ./output_momsolve --ls=1000. -t "ismipc_1000.toml"
 cp output_momsolve/ismipc_U_obs.h5 input/
 
 
@@ -18,6 +18,6 @@ RUN_DIR=$FENICS_ICE_BASE_DIR/runs/
 python $RUN_DIR/run_inv.py ismipc_1000_corr.toml
 python $RUN_DIR/run_forward.py ismipc_1000_corr.toml
 python $RUN_DIR/run_eigendec.py ismipc_1000_corr.toml
-#python $RUN_DIR/run_errorprop.py ismipc_1000_corr.toml
+python $RUN_DIR/run_errorprop.py ismipc_1000_corr.toml
 python $RUN_DIR/run_invsigma.py ismipc_1000_corr.toml
-#python $RUN_DIR/run_sample_post.py ismipc_1000_corr.toml
+python $RUN_DIR/run_sample_post.py ismipc_1000_corr.toml
diff --git a/example_cases/ismipc_2000/ismipc_2000.sh b/example_cases/ismipc_2000/ismipc_2000.sh
index 099afb7..1187bab 100755
--- a/example_cases/ismipc_2000/ismipc_2000.sh
+++ b/example_cases/ismipc_2000/ismipc_2000.sh
@@ -18,4 +18,4 @@ python $RUN_DIR/run_forward.py ismipc_2000.toml
 python $RUN_DIR/run_eigendec.py ismipc_2000.toml
 python $RUN_DIR/run_errorprop.py ismipc_2000.toml
 python $RUN_DIR/run_invsigma.py ismipc_2000.toml
-#python $RUN_DIR/run_sample_post.py ismipc_2000.toml
+python $RUN_DIR/run_sample_post.py ismipc_2000.toml
diff --git a/example_cases/ismipc_2000_corr/ismipc_2000_corr.sh b/example_cases/ismipc_2000_corr/ismipc_2000_corr.sh
index 718bebf..9ebb9b6 100755
--- a/example_cases/ismipc_2000_corr/ismipc_2000_corr.sh
+++ b/example_cases/ismipc_2000_corr/ismipc_2000_corr.sh
@@ -6,9 +6,9 @@ python $FENICS_ICE_BASE_DIR/aux/gen_rect_mesh.py -o ./input/momsolve_mesh.xml -x
 python $FENICS_ICE_BASE_DIR/aux/gen_rect_mesh.py -o ./input/ismip_mesh.xml -xmax 40000 -ymax 40000 -nx 30 -ny 30
 python $FENICS_ICE_BASE_DIR/aux/gen_ismipC_domain.py -o ./input/ismipc_input.h5 -L 40000 -nx 100 -ny 100 --reflect
 python $FENICS_ICE_BASE_DIR/runs/run_momsolve.py momsolve.toml
-python $FENICS_ICE_BASE_DIR/aux/Uobs_from_momsolve.py -i "U.h5" -o "ismipc_U_obs_nocorr.h5" -L 40000 -d ./output_momsolve --ls=2000. 
-cp output_momsolve/ismipc_U_obs_nocorr.h5 input/
-python $FENICS_ICE_BASE_DIR/aux/Uobs_from_momsolve_cov.py -i "U.h5" -o "ismipc_U_obs.h5" -L 40000 -d ./output_momsolve --ls=2000. -t "ismipc_2000.toml"
+#python $FENICS_ICE_BASE_DIR/aux/Uobs_from_momsolve.py -i "U.h5" -o "ismipc_U_obs_nocorr.h5" -L 40000 -d ./output_momsolve --ls=2000. 
+#cp output_momsolve/ismipc_U_obs_nocorr.h5 input/
+python $FENICS_ICE_BASE_DIR/aux/Uobs_from_momsolve.py -i "U.h5" -o "ismipc_U_obs.h5" -L 40000 -d ./output_momsolve --ls=2000. -t "ismipc_2000.toml"
 cp output_momsolve/ismipc_U_obs.h5 input/
 
 
diff --git a/example_cases/ismipc_4000/ismipc_4000.sh b/example_cases/ismipc_4000/ismipc_4000.sh
index 08a604d..3e47b17 100755
--- a/example_cases/ismipc_4000/ismipc_4000.sh
+++ b/example_cases/ismipc_4000/ismipc_4000.sh
@@ -18,4 +18,4 @@ python $RUN_DIR/run_forward.py ismipc_4000.toml
 python $RUN_DIR/run_eigendec.py ismipc_4000.toml
 python $RUN_DIR/run_errorprop.py ismipc_4000.toml
 python $RUN_DIR/run_invsigma.py ismipc_4000.toml
-#python $RUN_DIR/run_sample_post.py ismipc_4000.toml
+python $RUN_DIR/run_sample_post.py ismipc_4000.toml
diff --git a/example_cases/ismipc_4000_corr/ismipc_4000_corr.sh b/example_cases/ismipc_4000_corr/ismipc_4000_corr.sh
index bd7a356..ae7a178 100755
--- a/example_cases/ismipc_4000_corr/ismipc_4000_corr.sh
+++ b/example_cases/ismipc_4000_corr/ismipc_4000_corr.sh
@@ -6,9 +6,9 @@ python $FENICS_ICE_BASE_DIR/aux/gen_rect_mesh.py -o ./input/momsolve_mesh.xml -x
 python $FENICS_ICE_BASE_DIR/aux/gen_rect_mesh.py -o ./input/ismip_mesh.xml -xmax 40000 -ymax 40000 -nx 30 -ny 30
 python $FENICS_ICE_BASE_DIR/aux/gen_ismipC_domain.py -o ./input/ismipc_input.h5 -L 40000 -nx 100 -ny 100 --reflect
 python $FENICS_ICE_BASE_DIR/runs/run_momsolve.py momsolve.toml
-python $FENICS_ICE_BASE_DIR/aux/Uobs_from_momsolve.py -i "U.h5" -o "ismipc_U_obs_nocorr.h5" -L 40000 -d ./output_momsolve --ls=4000. 
-cp output_momsolve/ismipc_U_obs_nocorr.h5 input/
-python $FENICS_ICE_BASE_DIR/aux/Uobs_from_momsolve_cov.py -i "U.h5" -o "ismipc_U_obs.h5" -L 40000 -d ./output_momsolve --ls=4000. -t "ismipc_4000.toml"
+#python $FENICS_ICE_BASE_DIR/aux/Uobs_from_momsolve.py -i "U.h5" -o "ismipc_U_obs_nocorr.h5" -L 40000 -d ./output_momsolve --ls=4000. 
+#cp output_momsolve/ismipc_U_obs_nocorr.h5 input/
+python $FENICS_ICE_BASE_DIR/aux/Uobs_from_momsolve.py -i "U.h5" -o "ismipc_U_obs.h5" -L 40000 -d ./output_momsolve --ls=4000. -t "ismipc_4000.toml"
 cp output_momsolve/ismipc_U_obs.h5 input/
 
 
@@ -18,6 +18,6 @@ RUN_DIR=$FENICS_ICE_BASE_DIR/runs/
 python $RUN_DIR/run_inv.py ismipc_4000_corr.toml
 python $RUN_DIR/run_forward.py ismipc_4000_corr.toml
 python $RUN_DIR/run_eigendec.py ismipc_4000_corr.toml
-#python $RUN_DIR/run_errorprop.py ismipc_4000_corr.toml
+python $RUN_DIR/run_errorprop.py ismipc_4000_corr.toml
 python $RUN_DIR/run_invsigma.py ismipc_4000_corr.toml
-#python $RUN_DIR/run_sample_post.py ismipc_4000_corr.toml
+python $RUN_DIR/run_sample_post.py ismipc_4000_corr.toml
diff --git a/example_cases/ismipc_500/ismipc_500.sh b/example_cases/ismipc_500/ismipc_500.sh
index a505290..25b35f3 100755
--- a/example_cases/ismipc_500/ismipc_500.sh
+++ b/example_cases/ismipc_500/ismipc_500.sh
@@ -18,4 +18,4 @@ python $RUN_DIR/run_forward.py ismipc_500.toml
 python $RUN_DIR/run_eigendec.py ismipc_500.toml
 python $RUN_DIR/run_errorprop.py ismipc_500.toml
 python $RUN_DIR/run_invsigma.py ismipc_500.toml
-#python $RUN_DIR/run_sample_post.py ismipc_500.toml
+python $RUN_DIR/run_sample_post.py ismipc_500.toml
diff --git a/example_cases/ismipc_500_corr/ismipc_500_corr.sh b/example_cases/ismipc_500_corr/ismipc_500_corr.sh
index 35b5bc3..945ef16 100755
--- a/example_cases/ismipc_500_corr/ismipc_500_corr.sh
+++ b/example_cases/ismipc_500_corr/ismipc_500_corr.sh
@@ -6,9 +6,9 @@ python $FENICS_ICE_BASE_DIR/aux/gen_rect_mesh.py -o ./input/momsolve_mesh.xml -x
 python $FENICS_ICE_BASE_DIR/aux/gen_rect_mesh.py -o ./input/ismip_mesh.xml -xmax 40000 -ymax 40000 -nx 30 -ny 30
 python $FENICS_ICE_BASE_DIR/aux/gen_ismipC_domain.py -o ./input/ismipc_input.h5 -L 40000 -nx 100 -ny 100 --reflect
 python $FENICS_ICE_BASE_DIR/runs/run_momsolve.py momsolve.toml
-python $FENICS_ICE_BASE_DIR/aux/Uobs_from_momsolve.py -i "U.h5" -o "ismipc_U_obs_nocorr.h5" -L 40000 -d ./output_momsolve --ls=500. 
-cp output_momsolve/ismipc_U_obs_nocorr.h5 input/
-python $FENICS_ICE_BASE_DIR/aux/Uobs_from_momsolve_cov.py -i "U.h5" -o "ismipc_U_obs.h5" -L 40000 -d ./output_momsolve --ls=500. -t "ismipc_500.toml"
+#python $FENICS_ICE_BASE_DIR/aux/Uobs_from_momsolve.py -i "U.h5" -o "ismipc_U_obs_nocorr.h5" -L 40000 -d ./output_momsolve --ls=500. 
+#cp output_momsolve/ismipc_U_obs_nocorr.h5 input/
+python $FENICS_ICE_BASE_DIR/aux/Uobs_from_momsolve.py -i "U.h5" -o "ismipc_U_obs.h5" -L 40000 -d ./output_momsolve --ls=500. -t "ismipc_500.toml"
 cp output_momsolve/ismipc_U_obs.h5 input/
 
 
@@ -18,6 +18,6 @@ RUN_DIR=$FENICS_ICE_BASE_DIR/runs/
 python $RUN_DIR/run_inv.py ismipc_500_corr.toml
 python $RUN_DIR/run_forward.py ismipc_500_corr.toml
 python $RUN_DIR/run_eigendec.py ismipc_500_corr.toml
-#python $RUN_DIR/run_errorprop.py ismipc_500_corr.toml
+python $RUN_DIR/run_errorprop.py ismipc_500_corr.toml
 python $RUN_DIR/run_invsigma.py ismipc_500_corr.toml
-#python $RUN_DIR/run_sample_post.py ismipc_500_corr.toml
+python $RUN_DIR/run_sample_post.py ismipc_500_corr.toml
diff --git a/example_cases/ismipc_8000/ismipc_8000.sh b/example_cases/ismipc_8000/ismipc_8000.sh
index 2bfc4a7..01b2be2 100755
--- a/example_cases/ismipc_8000/ismipc_8000.sh
+++ b/example_cases/ismipc_8000/ismipc_8000.sh
@@ -18,4 +18,4 @@ python $RUN_DIR/run_forward.py ismipc_8000.toml
 python $RUN_DIR/run_eigendec.py ismipc_8000.toml
 python $RUN_DIR/run_errorprop.py ismipc_8000.toml
 python $RUN_DIR/run_invsigma.py ismipc_8000.toml
-#python $RUN_DIR/run_sample_post.py ismipc_8000.toml
+python $RUN_DIR/run_sample_post.py ismipc_8000.toml
diff --git a/example_cases/ismipc_8000_corr/ismipc_8000_corr.sh b/example_cases/ismipc_8000_corr/ismipc_8000_corr.sh
index e21168f..4d9f803 100755
--- a/example_cases/ismipc_8000_corr/ismipc_8000_corr.sh
+++ b/example_cases/ismipc_8000_corr/ismipc_8000_corr.sh
@@ -6,9 +6,9 @@ python $FENICS_ICE_BASE_DIR/aux/gen_rect_mesh.py -o ./input/momsolve_mesh.xml -x
 python $FENICS_ICE_BASE_DIR/aux/gen_rect_mesh.py -o ./input/ismip_mesh.xml -xmax 40000 -ymax 40000 -nx 30 -ny 30
 python $FENICS_ICE_BASE_DIR/aux/gen_ismipC_domain.py -o ./input/ismipc_input.h5 -L 40000 -nx 100 -ny 100 --reflect
 python $FENICS_ICE_BASE_DIR/runs/run_momsolve.py momsolve.toml
-python $FENICS_ICE_BASE_DIR/aux/Uobs_from_momsolve.py -i "U.h5" -o "ismipc_U_obs_nocorr.h5" -L 40000 -d ./output_momsolve --ls=8000. 
-cp output_momsolve/ismipc_U_obs_nocorr.h5 input/
-python $FENICS_ICE_BASE_DIR/aux/Uobs_from_momsolve_cov.py -i "U.h5" -o "ismipc_U_obs.h5" -L 40000 -d ./output_momsolve --ls=8000. -t "ismipc_8000.toml"
+#python $FENICS_ICE_BASE_DIR/aux/Uobs_from_momsolve.py -i "U.h5" -o "ismipc_U_obs_nocorr.h5" -L 40000 -d ./output_momsolve --ls=8000. 
+#cp output_momsolve/ismipc_U_obs_nocorr.h5 input/
+python $FENICS_ICE_BASE_DIR/aux/Uobs_from_momsolve.py -i "U.h5" -o "ismipc_U_obs.h5" -L 40000 -d ./output_momsolve --ls=8000. -t "ismipc_8000.toml"
 cp output_momsolve/ismipc_U_obs.h5 input/
 
 
@@ -18,6 +18,6 @@ RUN_DIR=$FENICS_ICE_BASE_DIR/runs/
 python $RUN_DIR/run_inv.py ismipc_8000_corr.toml
 python $RUN_DIR/run_forward.py ismipc_8000_corr.toml
 python $RUN_DIR/run_eigendec.py ismipc_8000_corr.toml
-#python $RUN_DIR/run_errorprop.py ismipc_8000_corr.toml
+python $RUN_DIR/run_errorprop.py ismipc_8000_corr.toml
 python $RUN_DIR/run_invsigma.py ismipc_8000_corr.toml
-#python $RUN_DIR/run_sample_post.py ismipc_8000_corr.toml
+python $RUN_DIR/run_sample_post.py ismipc_8000_corr.toml
diff --git a/example_cases/plots_gmd/make_plots.py b/example_cases/plots_gmd/make_plots.py
new file mode 100644
index 0000000..7d5f55b
--- /dev/null
+++ b/example_cases/plots_gmd/make_plots.py
@@ -0,0 +1,59 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import os
+
+# Lcurve
+
+dir='../ismipc_rc_1e6/output/'
+
+gam = np.load(dir + 'LCurvGam.npy')
+norm = np.load(dir + 'LCurvNorm.npy')
+mis = np.load(dir + 'LCurvMis.npy')
+
+fig,ax = plt.subplots()
+
+for i in range(len(gam)):
+  ax.plot(norm[i],mis[i],'o',markersize=4,markerfacecolor='blue',markeredgecolor='k')
+
+for i in range(len(gam)):
+  if (i==4):
+   ax.annotate(str(gam[i]),(norm[i]-6,mis[i]-6))
+  else:
+   ax.annotate(str(gam[i]),(norm[i]+1,mis[i]+1))
+
+plt.xlabel('Regularisation cost')
+plt.ylabel('Model-data misfit cost')
+
+
+plt.savefig('L_curve.png')
+
+plt.close(fig)
+
+pwd = os.getcwd()
+
+os.system('python ../../aux/plotting/plot_inv_results_hist.py ismipc_rc_1e6 ' + pwd)
+os.system('python ../../aux/plotting/plot_inv_results_hist.py ismipc_rc_1e4 ' + pwd)
+
+plt.close(fig)
+
+os.system('python ../../aux/plotting/plot_leading_eigenfuncs.py ' + pwd)
+
+plt.close(fig)
+
+os.system("python ../../aux/plotting/plotting_new/plot_eigenvalue_decay_gen_new.py $PWD reg_decay.png ismipc_rc_1e6 '$\gamma_{alpha}$=50' ismipc_rc_1e5 '$\gamma_{alpha}$=10' ismipc_rc_1e4 '$\gamma_{alpha}$=1'")
+
+plt.close(fig)
+
+os.system("python ../../aux/plotting/plotting_new/plot_eigenvalue_decay_gen_new.py $PWD res_decay.png ismipc_20x20 '$\Delta$ x = 2 km' ismipc_rc_1e5 '$\Delta$ x = 1.33 km' ismipc_40x40 '$\Delta$ x = 1 km'")
+
+plt.close(fig)
+
+os.system("python ../../aux/plotting/plotting_new/plot_eigenvalue_decay_gen_new.py $PWD sample_decay.png ismipc_8000 8000m ismipc_4000 4km ismipc_2000 2km ismipc_1000 1km ismipc_500 500m")
+
+plt.close(fig)
+
+#os.system("python ../../aux/plotting/plotting_new/plot_eigenvalue_decay_gen_new.py $PWD sample_decay_autocorr_alt2.png corr_alt/ismipc_8000_corr 8km corr_alt/ismipc_4000_corr 4km corr_alt/ismipc_2000_corr 2km corr_alt/ismipc_1000_corr 1km corr_alt/ismipc_500_corr 500m")
+
+#plt.close(fig)
+
+os.system('python ../../aux/plotting/plotting_new/plot_paths_sample.py $PWD paths.png')
diff --git a/example_cases/run_all.sh b/example_cases/run_all.sh
index 9c7f092..c300984 100644
--- a/example_cases/run_all.sh
+++ b/example_cases/run_all.sh
@@ -1,7 +1,7 @@
 #!/bin/bash
 #Runs each model in turn (after first clearing out any previous results)
 
-for f in $(ls -d ismipc*00 ismipc_rc*);
+for f in $(ls -d ismipc*corr);
 do echo $f
 cd $f
 echo "Running model in $f"
diff --git a/fenics_ice/model.py b/fenics_ice/model.py
index 4f197bd..f704283 100644
--- a/fenics_ice/model.py
+++ b/fenics_ice/model.py
@@ -28,6 +28,7 @@
 import logging
 from IPython import embed
 from scipy.linalg import inv as Invert
+from scipy.linalg import norm as matNorm
 
 log = logging.getLogger("fenics_ice")
 
@@ -77,6 +78,8 @@ def __init__(self, mesh_in, input_data, param_in, init_fields=True,
         self.def_lat_dirichletbc()
         self.GammaInvObsU = None
         self.GammaInvObsV = None
+        self.McovU = None
+        self.McovV = None
 
         if init_fields:
             self.init_fields_from_data()
@@ -262,7 +265,8 @@ def init_obs_inv_cov(self):
         xy = self.uv_obs_pts
         npts = np.size(xy,0)
         dist_decay = self.params.obs.pts_autocorr
-        Mcov = np.zeros((npts,npts))
+        self.McovU = np.zeros((npts,npts))
+        self.McovV = np.zeros((npts,npts))
         xpt1 = xy[:,0]
         ypt1 = xy[:,1]
         L = 40.e3
@@ -284,10 +288,13 @@ def init_obs_inv_cov(self):
           dist_func = np.sqrt(distx**2+disty**2)
 
           row = self.u_std[i] * self.u_std * np.exp(-dist_func/dist_decay)
-          Mcov[i,:] = row
+          self.McovU[i,:] = row
+
+#         mnorm = matNorm(self.McovU,2)
+#         self.McovU = self.McovU / mnorm
     
          t1 = time.perf_counter()
-         self.GammaInvObsU = Invert(Mcov)
+         self.GammaInvObsU = Invert(self.McovU)
          self.GammaInvObsU = .5*(self.GammaInvObsU + self.GammaInvObsU.T)
          t2 = time.perf_counter()
          print('U cov matrix invert, time= ' + str(t2-t1) + ' sec')
@@ -307,10 +314,13 @@ def init_obs_inv_cov(self):
           dist_func = np.sqrt(distx**2+disty**2)
 
           row = self.v_std[i] * self.v_std * np.exp(-dist_func/dist_decay)
-          Mcov[i,:] = row
+          self.McovV[i,:] = row
+
+#         mnorm = matNorm(self.McovV,2)
+#         self.McovV = self.McovV / mnorm
 
          t1 = time.perf_counter()
-         self.GammaInvObsV = Invert(Mcov)
+         self.GammaInvObsV = Invert(self.McovV)
          self.GammaInvObsV = .5*(self.GammaInvObsV + self.GammaInvObsV.T)
          t2 = time.perf_counter()
          print('V cov matrix invert, time= ' + str(t2-t1) + ' sec')
diff --git a/fenics_ice/prior.py b/fenics_ice/prior.py
index 4006502..bce6566 100644
--- a/fenics_ice/prior.py
+++ b/fenics_ice/prior.py
@@ -18,6 +18,7 @@
 from dolfin import *
 from tlm_adjoint import *
 from .decorators import count_calls, timer
+from IPython import embed
 from fenics_ice.sqrt_mass_matrix_action import A_root_action
 
 class laplacian(object):
diff --git a/fenics_ice/solver.py b/fenics_ice/solver.py
index 6c7ae56..a88a4de 100644
--- a/fenics_ice/solver.py
+++ b/fenics_ice/solver.py
@@ -772,16 +772,18 @@ def comp_J_inv(self, verbose=False, noMisfit=False, noReg=False):
          v_mismatch = ((v_pts-v_obs_pts)/v_std_pts)
          NormSqSolver(project(v_mismatch, obs_space), J_ls_term_v).solve()
         else:
-#         u_mismatch = Function(obs_space, name='u_mismatch')
+         u_mismatch = Function(obs_space, name='u_mismatch')
 #         u_mismatch.vector()[:] = u_pts.vector()[:]-u_obs_pts.vector()[:]
-         u_mismatch = project(u_pts-u_obs_pts,obs_space)
+         u_mis = u_pts-u_obs_pts
+         LocalProjectionSolver(u_mis,u_mismatch).solve()
          tmp_fcn_upts = Function(obs_space, name='tmp_fcn_upts')
          NumPyMatrixActionSolver(self.GammaInvObsU, u_mismatch, tmp_fcn_upts).solve()
          InnerProductSolver(u_mismatch, tmp_fcn_upts, J_ls_term_u).solve()
 
-#         v_mismatch = Function(obs_space, name='v_mismatch')
+         v_mismatch = Function(obs_space, name='v_mismatch')
 #         v_mismatch.vector()[:] = v_pts.vector()[:]-v_obs_pts.vector()[:]
-         v_mismatch = project(v_pts-v_obs_pts,obs_space)
+         v_mis = v_pts-v_obs_pts
+         LocalProjectionSolver(v_mis,v_mismatch).solve()
          tmp_fcn_vpts = Function(obs_space, name='tmp_fcn_vpts')
          NumPyMatrixActionSolver(self.GammaInvObsV, v_mismatch, tmp_fcn_vpts).solve()
          InnerProductSolver(v_mismatch, tmp_fcn_vpts, J_ls_term_v).solve()
diff --git a/fenics_ice/sqrt_mass_matrix_action.py b/fenics_ice/sqrt_mass_matrix_action.py
index cd8144d..d52c896 100644
--- a/fenics_ice/sqrt_mass_matrix_action.py
+++ b/fenics_ice/sqrt_mass_matrix_action.py
@@ -5,6 +5,7 @@
 
 import numpy as np
 import scipy.special as special
+from IPython import embed
 
 __all__ = \
     [
@@ -64,6 +65,7 @@ def A_action(x):
     y = x.copy()
     z = x.copy()
 
+
     j = 1
     while True:
         z.axpy(-1.0 / beta, A_action(z))
@@ -240,4 +242,4 @@ def test_sample_p1_mass():
     print(x_outer_mean)
     error_norm = abs(M_inv - x_outer_mean).max()
     print(f"error norm = {error_norm:.6e}")
-    assert error_norm < 0.09
\ No newline at end of file
+    assert error_norm < 0.09