diff --git a/example_cases/ismipc_30x30/ismipc_30x30.sh b/example_cases/ismipc_30x30/ismipc_30x30.sh
index fbbb1f1..6bc6ea3 100755
--- a/example_cases/ismipc_30x30/ismipc_30x30.sh
+++ b/example_cases/ismipc_30x30/ismipc_30x30.sh
@@ -13,8 +13,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_30x30.toml
-python $RUN_DIR/run_forward.py ismipc_30x30.toml
-python $RUN_DIR/run_eigendec.py ismipc_30x30.toml
-python $RUN_DIR/run_errorprop.py ismipc_30x30.toml
-python $RUN_DIR/run_invsigma.py ismipc_30x30.toml
+#python $RUN_DIR/run_inv.py ismipc_30x30.toml
+#python $RUN_DIR/run_forward.py ismipc_30x30.toml
+#python $RUN_DIR/run_eigendec.py ismipc_30x30.toml
+#python $RUN_DIR/run_errorprop.py ismipc_30x30.toml
+mpirun -n 4 python $RUN_DIR/run_obs_sens_prop.py ismipc_30x30.toml
diff --git a/example_cases/ismipc_30x30/ismipc_30x30.toml b/example_cases/ismipc_30x30/ismipc_30x30.toml
index 73afc83..60ce543 100644
--- a/example_cases/ismipc_30x30/ismipc_30x30.toml
+++ b/example_cases/ismipc_30x30/ismipc_30x30.toml
@@ -32,7 +32,7 @@ random_seed = 0
 [mesh]
 
 mesh_filename = "ismip_mesh.xml"
-periodic_bc = true
+periodic_bc = false
 
 [obs]
 
@@ -67,7 +67,8 @@ sliding_law = 'linear' #budd, linear
 [momsolve.picard_params]
 nonlinear_solver = "newton"
 [momsolve.picard_params.newton_solver]
-linear_solver = "umfpack"
+linear_solver = "cg"
+preconditioner = "hypre_amg"
 maximum_iterations = 200
 absolute_tolerance = 1.0e-0
 relative_tolerance = 1.0e-3
@@ -77,7 +78,9 @@ error_on_nonconvergence =  false
 [momsolve.newton_params]
 nonlinear_solver = "newton"
 [momsolve.newton_params.newton_solver]
-linear_solver = "umfpack"
+#linear_solver = "umfpack"
+linear_solver = "cg"
+preconditioner = "hypre_amg"
 maximum_iterations = 25
 absolute_tolerance = 1.0e-7
 relative_tolerance = 1.0e-8
@@ -128,6 +131,10 @@ qoi = 'h2' #or 'vaf'
     name = "Bottom Edge"
     id = 4
 
+[mass_solve]
+
+use_cg_thickness = true
+
 [testing]
 
 expected_init_alpha = 531.6114524861194
diff --git a/fenics_ice/config.py b/fenics_ice/config.py
index 6ff082e..680049b 100644
--- a/fenics_ice/config.py
+++ b/fenics_ice/config.py
@@ -124,6 +124,19 @@ def parse(self):
         except KeyError:
             pass
 
+        try:
+            obs_sens_dict = self.config_dict['obs_sens']
+        except KeyError:
+            obs_sens_dict = {}
+        self.obs_sens = ObsSensCfg(**obs_sens_dict)
+
+        try:
+            mass_solve_dict = self.config_dict['mass_solve']
+        except KeyError:
+            mass_solve_dict = {}
+        self.mass_solve = MassSolveCfg(**mass_solve_dict)
+
+
     def check_dirs(self):
         """
         Check input directory exists & create output dir if necessary.
@@ -138,13 +151,15 @@ def check_dirs(self):
                     self.time.phase_name,
                     self.eigendec.phase_name,
                     self.error_prop.phase_name,
-                    self.inv_sigma.phase_name]
+                    self.inv_sigma.phase_name,
+                    self.obs_sens.phase_name]
 
         ph_suffix = [self.inversion.phase_suffix,
                     self.time.phase_suffix,
                     self.eigendec.phase_suffix,
                     self.error_prop.phase_suffix,
-                    self.inv_sigma.phase_suffix]
+                    self.inv_sigma.phase_suffix,
+                    self.obs_sens.phase_suffix]
 
         for ph, suff in zip(ph_names, ph_suffix):
             out_dir = (outdir / ph / suff)
@@ -253,6 +268,17 @@ class ErrorPropCfg(ConfigPrinter):
     phase_name: str = 'error_prop'
     phase_suffix: str = ''
 
+
+@dataclass(frozen=True)
+class ObsSensCfg(ConfigPrinter):
+    """
+    Configuration related to observation sensitivities
+    """
+    qoi: str = 'vaf'
+    phase_name: str = 'obs_sens'
+    phase_suffix: str = ''
+    
+
 @dataclass(frozen=True)
 class SampleCfg(ConfigPrinter):
     """
@@ -394,6 +420,18 @@ def __post_init__(self):
             assert self.min_thickness >= 0.0
 
 
+@dataclass(frozen=True)
+class MassSolveCfg(ConfigPrinter):
+    """
+    Options for mass balance solver
+    """
+
+    use_cg_thickness: bool = False
+
+    def __post_init__(self):
+     """  """
+
+
 @dataclass(frozen=True)
 class MomsolveCfg(ConfigPrinter):
     """
@@ -513,7 +551,6 @@ def __post_init__(self):
 
         for fname in fname_default_suff:
             self.set_default_filename(fname, fname_default_suff[fname])
-            #embed()
 
 @dataclass(frozen=True)
 class TimeCfg(ConfigPrinter):
diff --git a/fenics_ice/inout.py b/fenics_ice/inout.py
index e2dc462..03ed8b2 100644
--- a/fenics_ice/inout.py
+++ b/fenics_ice/inout.py
@@ -294,7 +294,6 @@ def write_variable(var, params, name=None, outdir=None, phase_name='', phase_suf
     outvar.rename(name, "")
     # Prefix the run name
     outfname = Path(outdir) / phase_name / phase_suffix / "_".join((params.io.run_name+phase_suffix, name))
-    #embed()
 
     # Write out output according to user specified format in toml
     output_var_format = params.io.output_var_format
diff --git a/fenics_ice/model.py b/fenics_ice/model.py
index 064c164..d57e0cb 100644
--- a/fenics_ice/model.py
+++ b/fenics_ice/model.py
@@ -26,10 +26,43 @@
 from pathlib import Path
 from numpy.random import randn
 import logging
-from IPython import embed
 
 log = logging.getLogger("fenics_ice")
 
+        # Functions for repeated ungridded interpolation
+        # TODO - this will not handle extrapolation/missing data
+        # nicely - unfound simplex are returned '-1' which takes the last
+        # tri.simplices...
+
+        # at the moment i have moved these from vel_obs_from_data, so they 
+        # can be called directly from a run script.
+        # the ismipc test, which calls this function, still seems to perform fine
+        # but this refactoring may make things less efficient.
+def interp_weights(xy, uv, periodic_bc, d=2):
+            """Compute the nearest vertices & weights (for reuse)"""
+            from scipy.spatial import Delaunay
+            tri = Delaunay(xy)
+            simplex = tri.find_simplex(uv)
+
+            if not np.all(simplex >= 0):
+                if not periodic_bc:
+                    log.warning("Some points missing in interpolation "
+                              "of velocity obs to function space.")
+                else:
+                    log.warning("Some points missing in interpolation "
+                                "of velocity obs to function space.")
+
+            vertices = np.take(tri.simplices, simplex, axis=0)
+            temp = np.take(tri.transform, simplex, axis=0)
+            delta = uv - temp[:, d]
+            bary = np.einsum('njk,nk->nj', temp[:, :d, :], delta)
+            return vertices, np.hstack((bary, 1 - bary.sum(axis=1,
+                                                           keepdims=True)))
+
+def interpolate(values, vtx, wts):
+            """Bilinear interpolation, given vertices & weights above"""
+            return np.einsum('nj,nj->n', np.take(values, vtx), wts)
+
 class model:
     """
     The 'model' object is the core of any fenics_ice simulation. It handles loading input
@@ -122,7 +155,10 @@ def init_fields_from_data(self):
         self.bed = self.field_from_data("bed", self.Q)
         self.bmelt = self.field_from_data("bmelt", self.M, default=0.0)
         self.smb = self.field_from_data("smb", self.M, default=0.0)
-        self.H_np = self.field_from_data("thick", self.M, min_val=min_thick)
+        if (self.params.mass_solve.use_cg_thickness and self.params.mesh.periodic_bc):
+         self.H_np = self.field_from_data("thick", self.Qp, min_val=min_thick)
+        else:
+         self.H_np = self.field_from_data("thick", self.M, min_val=min_thick)
 
         if self.params.melt.use_melt_parameterisation:
        
@@ -308,34 +344,6 @@ def vel_obs_from_data(self):
                                use_cloud_point=self.params.inversion.use_cloud_point_velocities)
         else:
             inout.read_vel_obs(infile, model=self)
-        # Functions for repeated ungridded interpolation
-        # TODO - this will not handle extrapolation/missing data
-        # nicely - unfound simplex are returned '-1' which takes the last
-        # tri.simplices...
-        def interp_weights(xy, uv, d=2):
-            """Compute the nearest vertices & weights (for reuse)"""
-            from scipy.spatial import Delaunay
-            tri = Delaunay(xy)
-            simplex = tri.find_simplex(uv)
-
-            if not np.all(simplex >= 0):
-                if not self.params.mesh.periodic_bc:
-                    log.warning("Some points missing in interpolation "
-                              "of velocity obs to function space.")
-                else:
-                    log.warning("Some points missing in interpolation "
-                                "of velocity obs to function space.")
-
-            vertices = np.take(tri.simplices, simplex, axis=0)
-            temp = np.take(tri.transform, simplex, axis=0)
-            delta = uv - temp[:, d]
-            bary = np.einsum('njk,nk->nj', temp[:, :d, :], delta)
-            return vertices, np.hstack((bary, 1 - bary.sum(axis=1,
-                                                           keepdims=True)))
-
-        def interpolate(values, vtx, wts):
-            """Bilinear interpolation, given vertices & weights above"""
-            return np.einsum('nj,nj->n', np.take(values, vtx), wts)
 
         # Grab coordinates of both Lagrangian & DG function spaces
         # and compute (once) the interpolating arrays
@@ -343,9 +351,9 @@ def interpolate(values, vtx, wts):
         M_coords = self.M.tabulate_dof_coordinates()
 
         vtx_Q, wts_Q = interp_weights(self.vel_obs['uv_comp_pts'],
-                                      Q_coords)
+                                      Q_coords, self.params.mesh.periodic_bc)
         vtx_M, wts_M = interp_weights(self.vel_obs['uv_comp_pts'],
-                                      M_coords)
+                                      M_coords, self.params.mesh.periodic_bc)
 
         # Define new functions to hold results
         self.u_obs_Q = Function(self.Q, name="u_obs")
@@ -378,7 +386,7 @@ def interpolate(values, vtx, wts):
         # We need to do the same as above but for cloud point data
         # so we can write out a nicer output in the mesh coordinates
         vtx_Q_c, wts_Q_c = interp_weights(self.vel_obs['uv_obs_pts'],
-                                      Q_coords)
+                                      Q_coords, self.params.mesh.periodic_bc)
 
         # Define new functions to hold results
         self.u_cloud_Q = Function(self.Q, name="u_obs_cloud")
diff --git a/fenics_ice/solver.py b/fenics_ice/solver.py
index 5f93527..7c9a6cc 100644
--- a/fenics_ice/solver.py
+++ b/fenics_ice/solver.py
@@ -29,7 +29,6 @@
 import time
 import ufl
 import weakref
-from IPython import embed
 
 log = logging.getLogger("fenics_ice")
 
@@ -77,11 +76,26 @@ def interpolation_matrix(x_coords, y_space):
     return x_local, P
 
 
+def Amat_obs_action(P, Rvec, vec_cg, dg_space):
+    # This function implements the Rvec*P*D action on a P1 function
+    #  where D is a projection into DG space
+    # 
+    # to be called for each component of velocity
+    #
+
+    test, trial = TestFunction(dg_space), TrialFunction(dg_space)
+    vec_dg = Function(dg_space)
+    solve(inner(trial, test) * dx == inner(vec_cg, test) * dx,
+          vec_dg, solver_parameters={"linear_solver": "lu"})
+
+    return Rvec * (P @ vec_dg.vector().get_local())
+    
+
 class ssa_solver:
     """
     The ssa_solver object is currently the only kind of fenics_ice 'solver' available.
     """
-    def __init__(self, model, mixed_space=False):
+    def __init__(self, model, mixed_space=False, obs_sensitivity=False):
 
         # Enable aggressive compiler options
         parameters["form_compiler"]["optimize"] = False
@@ -93,6 +107,7 @@ def __init__(self, model, mixed_space=False):
         self.model.solvers.append(self)
         self.params = model.params
         self.mixed_space = mixed_space
+        self.obs_sensitivity = obs_sensitivity
 
         # Mesh/Function Spaces
         self.mesh = model.mesh
@@ -146,10 +161,15 @@ def __init__(self, model, mixed_space=False):
         self.U = Function(self.V, name="U")
         self.U_np = Function(self.V, name="U_np")
         self.Phi = TestFunction(self.V)
-        self.Ksi = TestFunction(self.M)
         self.pTau = TestFunction(self.Qp)
 
-        self.trial_H = TrialFunction(self.M)
+        if not (self.params.mass_solve.use_cg_thickness and self.params.mesh.periodic_bc):
+         self.trial_H = TrialFunction(self.M)
+         self.Ksi = TestFunction(self.M)
+        else:
+         self.trial_H = TrialFunction(self.Qp)
+         self.Ksi = TestFunction(self.Qp)
+
 
         # Facets
         self.ff = model.ff
@@ -607,14 +627,6 @@ def def_thickadv_eq(self):
             + inner(jump(Ksi), jump(0.5 * (dot(U_np, nm) + abs(dot(U_np, nm))) * trial_H))
             * dS
 
-            # Outflow at boundaries
-            + conditional(dot(U_np, nm) > 0, 1.0, 0.0)*inner(Ksi, dot(U_np * trial_H, nm))
-            * ds
-
-            # Inflow at boundaries
-            + conditional(dot(U_np, nm) < 0, 1.0, 0.0)*inner(Ksi, dot(U_np * H_init, nm))
-            * ds
-
             # basal melting
             + bmelt*Ksi*dx
 
@@ -622,6 +634,19 @@ def def_thickadv_eq(self):
             - smb*Ksi*dx
         )
 
+
+        if not (self.params.mass_solve.use_cg_thickness and self.params.mesh.periodic_bc):
+            self.thickadv = self.thickadv + (
+
+            # Outflow at boundaries
+                + conditional(dot(U_np, nm) > 0, 1.0, 0.0)*inner(Ksi, dot(U_np * trial_H, nm))
+                * ds
+
+           # Inflow at boundaries
+                + conditional(dot(U_np, nm) < 0, 1.0, 0.0)*inner(Ksi, dot(U_np * H_init, nm))
+                * ds
+            )
+
         # # Forward euler
         # self.thickadv = (inner(Ksi, ((trial_H - H_np) / dt)) * dx
         # - inner(grad(Ksi), U_np * H_np) * dx
@@ -1303,9 +1328,16 @@ def comp_J_inv(self, verbose=False):
                 J_v_obs, op=MPI.SUM)
             J_v_obs, = J_v_obs
 
+            u_std_local = u_std[obs_local]
+            v_std_local = v_std[obs_local]
+
             self._cached_J_mismatch_data \
                 = (interp_space,
                    u_PRP, v_PRP, l_u_obs, l_v_obs, J_u_obs, J_v_obs)
+            if (self.obs_sensitivity):
+                self._cached_Amat_vars = \
+                    (P, u_std_local, v_std_local, obs_local, interp_space)
+
         (interp_space,
          u_PRP, v_PRP, l_u_obs, l_v_obs, J_u_obs, J_v_obs) = \
             self._cached_J_mismatch_data
diff --git a/runs/run_forward.py b/runs/run_forward.py
index baddfd3..1d4c8de 100644
--- a/runs/run_forward.py
+++ b/runs/run_forward.py
@@ -72,6 +72,7 @@ def run_forward(config_file):
 
     # Run the forward model
     Q = slvr.timestep(adjoint_flag=1, qoi_func=qoi_func)
+
     # Run the adjoint model, computing gradient of Qoi w.r.t cntrl
     dQ_ts = compute_gradient(Q, cntrl)  # Isaac 27
 
diff --git a/runs/run_obs_sens_prop.py b/runs/run_obs_sens_prop.py
new file mode 100644
index 0000000..7ececf6
--- /dev/null
+++ b/runs/run_obs_sens_prop.py
@@ -0,0 +1,278 @@
+# For fenics_ice copyright information see ACKNOWLEDGEMENTS in the fenics_ice
+# root directory
+
+# This file is part of fenics_ice.
+#
+# fenics_ice is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published by
+# the Free Software Foundation, version 3 of the License.
+#
+# fenics_ice is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with tlm_adjoint.  If not, see <https://www.gnu.org/licenses/>.
+
+from fenics_ice.backend import Function, HDF5File
+from tlm_adjoint import set_manager, stop_manager, \
+        configure_tlm, function_tlm, restore_manager, \
+        EquationManager, start_manager
+
+import os
+
+import mpi4py.MPI as MPI  # noqa: N817
+from pathlib import Path
+import pickle
+import numpy as np
+import sys
+
+from fenics_ice import model, solver, inout
+from fenics_ice import mesh as fice_mesh
+from fenics_ice.config import ConfigParser
+from ufl import split
+from fenics_ice.solver import Amat_obs_action
+
+import matplotlib as mpl
+mpl.use("Agg")
+from scipy.sparse import spdiags
+
+os.environ["OMP_NUM_THREADS"] = "1"
+os.environ["OPENBLAS_NUM_THREADS"] = "1"
+
+
+# the function below is not in run_errorprop.py, but needed for obs sens
+@restore_manager
+def compute_tau(forward, u, m, dm):
+    # this block of code will do the "forward" (calculation of velocity and cost function) once
+    # and then find the jacobian of u in the direction needed
+    set_manager(EquationManager(cp_method="none", cp_parameters={}))
+    stop_manager()
+
+    start_manager(tlm=True)
+    configure_tlm((m, dm))
+    forward(m)
+    return function_tlm(u, (m, dm))        
+
+
+def run_obs_sens_prop(config_file):
+
+    # Read run config file
+    params = ConfigParser(config_file)
+    num_sens = params.time.num_sens
+
+    # Setup logging
+    inout.setup_logging(params)
+    inout.log_preamble("obssensprop", params)
+
+    outdir = params.io.output_dir
+
+    # Load the static model data (geometry, smb, etc)
+    input_data = inout.InputData(params)
+
+    # Eigen decomposition params
+    phase_eigen = params.eigendec.phase_name
+    phase_suffix_e = params.eigendec.phase_suffix
+    lamfile = params.io.eigenvalue_file
+    vecfile = params.io.eigenvecs_file
+
+    # Qoi forward params
+    phase_time = params.time.phase_name
+    phase_suffix_qoi = params.time.phase_suffix
+    dqoi_h5file = params.io.dqoi_h5file
+
+    if len(phase_suffix_e) > 0:
+        lamfile = params.io.run_name + phase_suffix_e + '_eigvals.p'
+        vecfile = params.io.run_name + phase_suffix_e + '_vr.h5'
+    if len(phase_suffix_qoi) > 0:
+        dqoi_h5file = params.io.run_name + phase_suffix_qoi + '_dQ_ts.h5'
+
+    # Get model mesh
+    mesh = fice_mesh.get_mesh(params)
+
+    # Define the model
+    mdl = model.model(mesh, input_data, params)
+
+    # Load alpha/beta fields
+    mdl.alpha_from_inversion()
+    mdl.beta_from_inversion()
+    mdl.bglen_from_data(mask_only=True)
+
+    # Setup our solver object -- obs_sensitivity flag set to ensure caching
+    slvr = solver.ssa_solver(mdl, mixed_space=params.inversion.dual, obs_sensitivity=True)
+
+    cntrl = slvr.get_control()
+    space = slvr.get_control_space()
+        
+    # call to slvr.forward() needed to cache variables
+    slvr.forward(cntrl)
+
+    # Regularization operator using inversion delta/gamma values
+    Prior = mdl.get_prior()
+    reg_op = Prior(slvr, space)
+
+    # Loads eigenvalues from file
+    outdir_e = Path(outdir)/phase_eigen/phase_suffix_e
+    with open(outdir_e/lamfile, 'rb') as ff:
+        eigendata = pickle.load(ff)
+        lam = eigendata[0].real.astype(np.float64)
+        nlam = len(lam)
+
+    # Check if eigendecomposition successfully produced num_eig
+    # or if some are NaN
+    if np.any(np.isnan(lam)):
+        raise RuntimeError("NaN eigenvalue(s)")
+
+    # and eigenvectors from .h5 file
+    eps = params.constants.float_eps
+    W = []
+    with HDF5File(MPI.COMM_WORLD, str(outdir_e/vecfile), 'r') as hdf5data:
+        for i in range(nlam):
+            w = Function(space)
+            hdf5data.read(w, f'v/vector_{i}')
+
+            # Test squared norm in prior == 1.0
+            B_inv_w = Function(space, space_type="conjugate_dual")
+            reg_op.action(w.vector(), B_inv_w.vector())
+            norm_sq_in_prior = w.vector().inner(B_inv_w.vector())
+            assert (abs(norm_sq_in_prior - 1.0) < eps)
+            del B_inv_w
+
+            W.append(w)
+
+    D = np.diag(lam / (lam + 1))  # D_r Isaac 20
+
+    # File containing dQoi_dCntrl (i.e. Jacobian of parameter to observable (Qoi))
+    outdir_qoi = Path(outdir)/phase_time/phase_suffix_qoi
+    hdf5data = HDF5File(MPI.COMM_WORLD, str(outdir_qoi/dqoi_h5file), 'r')
+
+    dQ_cntrl = Function(space, space_type="conjugate_dual")
+
+    run_length = params.time.run_length
+    num_sens = params.time.num_sens
+    t_sens = np.flip(np.linspace(run_length, 0, num_sens))
+
+    # above this point there is very little difference with run_errorprop.py
+    # -- exceptions are the defined function compute_tau() and 
+    # -- the call to solver.forward() to generate interpolation matrix
+
+    # below this point the result
+    #
+    # dQ/d\hat{p} = A (dp/dm) H^{-1} dQ_T/dm
+    #  
+    # is calculated, where p and \hat{p} are components of velocity (modelled and observed, respectively)
+    # Q_T is the QoI at a certain time, and A is a product of several operators involving interpolation, 
+    # projection, and observational covariance.
+
+    # this simply loads cached variables from solver.py which are cached in 
+    #  the above forward() call
+    assert hasattr(slvr,'_cached_Amat_vars'),\
+        "Attempt to retrieve Amat from solver type without first caching"
+    (P, u_std_local, v_std_local, obs_local, interp_space) = \
+                slvr._cached_Amat_vars
+
+    # This is the matrix R in the definition of matrix A (separate for each u and v)
+    Ru = spdiags(1.0 / (u_std_local ** 2),
+                              0, P.shape[0], P.shape[0])
+    Rv = spdiags(1.0 / (v_std_local ** 2),
+                              0, P.shape[0], P.shape[0])
+
+    dObsU = []
+    dObsV = []
+
+    comm = MPI.COMM_WORLD
+    size = comm.Get_size()
+    rank = comm.Get_rank()
+
+    for j in range(num_sens):
+
+        # for each time level T, we have a Q_T, hence the loop
+            
+        hdf5data.read(dQ_cntrl, f'dQd{cntrl[0].name()}/vector_{j}')
+
+        # the rest of this loop implements the same operations as 
+        # lines 137-147 of run_errorprop.py, ie
+        #
+        #   (Gamma_{prior} - W D W^T) acting on (dQ/dm)
+        #
+        # P3 needed to be defined as it is because
+        # P_vec = P2.vector() - P1.vector() led to errors later on (cannot recall specifics)
+            
+        tmp1 = np.asarray([w.vector().inner(dQ_cntrl.vector()) for w in W])
+        tmp2 = np.dot(D, tmp1)
+
+        P1 = Function(space)
+        for tmp, w in zip(tmp2, W):
+            P1.vector().axpy(tmp, w.vector())
+
+        P2 = Function(space)
+        reg_op.inv_action(dQ_cntrl.vector(), P2.vector())
+
+        P3 = Function(space)
+        P3.vector()[:] = P2.vector()[:] - P1.vector()[:]
+        
+
+        # tau is  d(U,V)/dm * (Gamma_{prior} - W D W^T) * (dQ/dm), or
+        #         d(U,V)/dm * P3
+        tau = compute_tau(slvr.forward, slvr.U, cntrl, P3)
+
+        # tau is in the space of U (right?)
+        tauu, tauv = split(tau)
+
+
+        # this block of code then implements A.tau
+        # but it needs to be made negative..
+
+        # note -- added mult by -1 before, but this was found to be incorrect
+        dobsu = Amat_obs_action(P, Ru, tauu, interp_space)
+        dobsv = Amat_obs_action(P, Rv, tauv, interp_space)
+
+        # this end result (above) corresponds only to the velocity obs
+        # that live on this processor's subdomain. An MPI reduce
+        # is used to generate the global vectors, which are then 
+        # saved via pickle
+
+        # the broadcast below is necessary to interpolate to a DG function
+        # but would not be needed o/w since the pickle dump is only on p0
+            
+        sendbuf = np.zeros(len(mdl.vel_obs['u_obs']),dtype='float')
+        rcvbuf = None
+        if rank == 0:
+            rcvbuf = np.empty(len(mdl.vel_obs['u_obs']),dtype='float')
+        sendbuf[obs_local] = dobsu
+        comm.Reduce(sendbuf, rcvbuf, root=0, op=MPI.SUM)
+        dobsU_0 = rcvbuf
+        dobsU = comm.bcast(dobsU_0, root=0)
+
+        sendbuf[:]=0.
+        if rank == 0:
+            rcvbuf = np.empty(len(mdl.vel_obs['u_obs']),dtype='float')
+        sendbuf[obs_local] = dobsv
+        comm.Reduce(sendbuf, rcvbuf, root=0, op=MPI.SUM)
+        dobsV_0 = rcvbuf
+        dobsV = comm.bcast(dobsV_0, root=0)
+
+        dObsU.append(dobsU)
+        dObsV.append(dobsV)
+
+
+    # Save data in diagnostics
+    phase_sens = params.obs_sens.phase_name
+    phase_suffix_sens = params.obs_sens.phase_suffix
+    diag_dir = Path(params.io.diagnostics_dir)/phase_sens/phase_suffix_sens
+
+    if rank == 0:
+        obs_pts = mdl.vel_obs['uv_obs_pts']
+        uobs = mdl.vel_obs['u_obs']
+        vobs = mdl.vel_obs['v_obs']
+
+        obs_sens_arrays = {'num_sens': num_sens, 'uv_obs_pts': obs_pts, \
+          'dObsU': dObsU, 'dObsV': dObsV, 'u_obs': uobs, 'v_obs': vobs}
+
+        with open(os.path.join(diag_dir, 'vel_sens.pkl'), 'wb') as f:
+          pickle.dump(obs_sens_arrays, f)
+ 
+if __name__ == "__main__":
+    assert len(sys.argv) == 2, "Expected a configuration file (*.toml)"
+    run_obs_sens_prop(sys.argv[1])