diff --git a/src/soca/Covariance/soca_covariance_mod.F90 b/src/soca/Covariance/soca_covariance_mod.F90
index c0446e741..3ed0c615c 100644
--- a/src/soca/Covariance/soca_covariance_mod.F90
+++ b/src/soca/Covariance/soca_covariance_mod.F90
@@ -1,4 +1,4 @@
-! (C) Copyright 2017-2021 UCAR.
+! (C) Copyright 2017-2023 UCAR.
 !
 ! This software is licensed under the terms of the Apache Licence Version 2.0
 ! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
@@ -34,6 +34,7 @@ module soca_covariance_mod
    real(kind=kind_real), allocatable :: pert_scale(:) !< index matches "vars"
    type(oops_variables), allocatable :: conv_vars(:)  !< index mathces "conv"
 
+   type(soca_geom), pointer :: geom
 contains
   !> \copybrief soca_cov_setup \see soca_cov_setup
   procedure :: setup => soca_cov_setup
@@ -67,7 +68,7 @@ module soca_covariance_mod
 subroutine soca_cov_setup(self, f_conf, geom, bkg, vars)
   class(soca_cov),        intent(inout) :: self   !< The covariance structure
   type(fckit_configuration), intent(in) :: f_conf !< The configuration
-  type(soca_geom),           intent(in) :: geom   !< Geometry
+  type(soca_geom), pointer,  intent(in) :: geom   !< Geometry
   type(soca_state),  target, intent(in) :: bkg    !< Background
   type(oops_variables),      intent(in) :: vars   !< List of variables
 
@@ -77,6 +78,9 @@ subroutine soca_cov_setup(self, f_conf, geom, bkg, vars)
   character(len=:), allocatable :: domain
   integer :: i, isc, iec, jsc, jec, ivar
 
+  ! save a copy of geometry for use later
+  self%geom => geom
+
   ! Setup list of variables to apply B on
   self%vars = vars
 
@@ -181,6 +185,9 @@ subroutine soca_cov_C_mult(self, dx)
   integer :: i, z
   type(soca_field), pointer :: field
   type(bump_type), pointer :: conv
+  real(kind=kind_real), allocatable :: fld_2d(:,:)
+
+  allocate(fld_2d(self%geom%isd:self%geom%ied, self%geom%jsd:self%geom%jed))
 
   do i = 1, self%vars%nvars()
     ! why is this sometimes getting an "empty" list with "none" in it?
@@ -196,7 +203,9 @@ subroutine soca_cov_C_mult(self, dx)
 
     ! apply convolution on each level
     do z = 1, field%nz
-      call soca_2d_convol(field%val(:,:,z), conv, dx%geom)
+      fld_2d = field%val(:,:,z)
+      call soca_2d_convol(fld_2d, conv, dx%geom)
+      field%val(:,:,z) = fld_2d      
     end do
   end do
 end subroutine soca_cov_C_mult
@@ -214,6 +223,9 @@ subroutine soca_cov_sqrt_C_mult(self, dx)
   type(soca_field), pointer :: field
   real(kind=kind_real) :: scale
   type(bump_type), pointer :: conv
+  real(kind=kind_real), allocatable :: fld_2d(:,:)
+
+  allocate(fld_2d(self%geom%isd:self%geom%ied, self%geom%jsd:self%geom%jed))
 
   do i = 1, self%vars%nvars()
     conv => null()
@@ -236,7 +248,9 @@ subroutine soca_cov_sqrt_C_mult(self, dx)
 
     ! apply convolution
     do z = 1,field%nz
-      call soca_2d_sqrt_convol(field%val(:,:,z), conv, dx%geom, scale)
+      fld_2d = field%val(:,:,z)
+      call soca_2d_sqrt_convol(fld_2d, conv, dx%geom, scale)
+      field%val(:,:,z) = fld_2d
     end do
 
   end do
@@ -264,49 +278,25 @@ subroutine soca_bump_correlation(self, horiz_convol, geom, f_conf_bump, f_conf_d
   integer :: i
   integer, allocatable :: hmask(:,:)
   integer, pointer :: int_ptr(:,:)
-  real(kind=kind_real), pointer :: real_ptr(:,:)
-  real(kind=kind_real), allocatable :: area(:)
+  real(kind=kind_real), pointer :: real_ptr(:,:), vArea(:,:), vRossby(:,:)
   type(atlas_functionspace) :: afunctionspace
   type(fieldset_type) :: afieldset, rh, rv
-  type(atlas_field) :: afield
+
+  type(atlas_field) :: afield, fArea, fRossby
   real(kind=kind_real) :: r_base, r_mult, r_min, r_max, r_min_grid
 
-  ! Wrap functionspace
-  afunctionspace = atlas_functionspace(geom%functionspaceInchalo%c_ptr())
+  ! wrap the functionspace
+  afunctionspace = atlas_functionspace(geom%functionspace%c_ptr())
 
   ! Geometry fieldset setup
   afieldset = atlas_fieldset()
 
-  ! Add area
-  afield = geom%functionspaceInchalo%create_field(name='area', kind=atlas_real(kind_real), levels=1)
-  call afield%data(real_ptr)
-  area = pack(geom%cell_area,geom%valid_halo_mask)
-  real_ptr(1,:) = area
-  call afieldset%add(afield)
-  call afield%final()
-
-  ! Add vertical unit
-  afield = geom%functionspaceInchalo%create_field(name='vert_coord', kind=atlas_real(kind_real), levels=1)
-  call afield%data(real_ptr)
-  real_ptr(1,:) = 1.0
-  call afieldset%add(afield)
-  call afield%final()
-
-  ! Add geographical mask
-  afield = geom%functionspaceInchalo%create_field(name='gmask', kind=atlas_integer(kind(0)), levels=1)
-  call afield%data(int_ptr)
-  int_ptr(1,:) = int(pack(geom%mask2d,geom%valid_halo_mask))
-  call afieldset%add(afield)
-  call afield%final()
-
-  afield = geom%functionspaceInchalo%create_field(name='owned', kind=atlas_integer(kind(0)), levels=1)
-  allocate(hmask(geom%isd:geom%ied, geom%jsd:geom%jed))
-  hmask = 0
-  hmask(geom%isc:geom%iec, geom%jsc:geom%jec) = 1
-  call afield%data(int_ptr)
-  int_ptr(1,:) = pack(hmask, geom%valid_halo_mask)
-  call afieldset%add(afield)
-  call afield%final()
+  ! add existing fields that were created by geometry
+  call afieldset%add(geom%fieldset%field('area'))
+  call afieldset%add(geom%fieldset%field('vert_coord'))
+  call afieldset%add(geom%fieldset%field('gmask'))
+  call afieldset%add(geom%fieldset%field('owned'))
+
 
   ! Set verbosity
   horiz_convol%mpl%verbose = (geom%f_comm%rank()==0)
@@ -315,6 +305,12 @@ subroutine soca_bump_correlation(self, horiz_convol, geom, f_conf_bump, f_conf_d
   call horiz_convol%create(geom%f_comm,afunctionspace,afieldset,f_conf_bump)
 
   if (horiz_convol%nam%new_nicas) then
+    ! get some fields from the geom fieldset
+    fArea = geom%fieldset%field('area')
+    call fArea%data(vArea)
+    fRossby = geom%fieldset%field('rossby_radius')
+    call fRossby%data(vRossby)
+
     ! get parameters for correlation lengths
     if (.not. f_conf_domain%get('base value', r_base))    r_base = 0.0
     if (.not. f_conf_domain%get('rossby mult', r_mult))   r_mult = 0.0
@@ -328,13 +324,13 @@ subroutine soca_bump_correlation(self, horiz_convol, geom, f_conf_bump, f_conf_d
     ! 3) min/max are imposed based on "min value" and "max value"
     ! 4) converted from a gaussian sigma to Gaspari-Cohn cutoff distance
     rh = atlas_fieldset()
-    afield = geom%functionspaceInchalo%create_field('var',kind=atlas_real(kind_real),levels=1)
+    afield = geom%functionspace%create_field('var',kind=atlas_real(kind_real),levels=1)
     call rh%add(afield)
     call afield%data(real_ptr)
-    real_ptr(1,:) = r_base + r_mult*pack(geom%rossby_radius, geom%valid_halo_mask)
+    real_ptr(1,:) = r_base + r_mult*vRossby(1,:)
     ! min based on grid size
     if (r_min_grid .gt. 0.0) then
-      real_ptr(1,:) = max(real_ptr(1,:),  sqrt(area)*r_min_grid )
+      real_ptr(1,:) = max(real_ptr(1,:),  sqrt(vArea(1,:))*r_min_grid )
     end if
     real_ptr(1,:) = min(r_max, real_ptr(1,:))
     real_ptr(1,:) = max(r_min, real_ptr(1,:))
@@ -344,7 +340,7 @@ subroutine soca_bump_correlation(self, horiz_convol, geom, f_conf_bump, f_conf_d
 
      ! rv
      rv = atlas_fieldset()
-     afield = geom%functionspaceInchalo%create_field('var',kind=atlas_real(kind_real),levels=1)
+     afield = geom%functionspace%create_field('var',kind=atlas_real(kind_real),levels=1)
      call rv%add(afield)
      call afield%data(real_ptr)
      real_ptr = 1.0
@@ -357,6 +353,8 @@ subroutine soca_bump_correlation(self, horiz_convol, geom, f_conf_bump, f_conf_d
      ! Clean up
      call rh%final()
      call rv%final()
+     call fRossby%final()
+     call fArea%final()
   end if
 
   ! Run BUMP drivers
@@ -371,24 +369,41 @@ end subroutine soca_bump_correlation
 !! Used by soca_cov::mult()
 !! \relates soca_covariance_mod::soca_cov
 subroutine soca_2d_convol(dx, horiz_convol, geom)
-  real(kind=kind_real), intent(inout) :: dx(:,:)
+  real(kind=kind_real), allocatable, intent(inout) :: dx(:,:)
   type(bump_type),      intent(inout) :: horiz_convol
   type(soca_geom),         intent(in) :: geom
 
-  type(fieldset_type) :: tmp_incr
-
-  ! Allocate ATLAS tmp_increment and make copy of dx
-  call geom%struct2atlas(dx(:,:), tmp_incr)
+  type(fieldset_type) :: fieldset
+  type(atlas_field)   :: field
+  real(kind_real), pointer :: real_ptr(:,:)
+  integer :: i, j, n
+
+  ! array to atlas
+  ! (Yeah, this code is duplicated in a few places, but this whole
+  ! class is going away "soon" so I don't care)
+  fieldset = atlas_fieldset()
+  field = geom%functionspace%create_field('var', kind=atlas_real(kind_real), levels=1)
+  call fieldset%add(field)
+  call field%data(real_ptr)
+  do j=geom%jsc,geom%jec
+    do i=geom%isc,geom%iec
+      real_ptr(1,geom%atlas_ij2idx(i,j)) = dx(i,j)
+    end do
+  end do
 
   ! Apply 2D convolution
-  call horiz_convol%apply_nicas(tmp_incr)
-
-  ! Copy ATLAS tmp_incr to structured dx
-  call geom%atlas2struct(dx(:,:), tmp_incr)
+  call horiz_convol%apply_nicas(fieldset)
 
+  ! atlas to array
+  do j=geom%jsc,geom%jec
+    do i=geom%isc,geom%iec
+      dx(i,j) = real_ptr(1,geom%atlas_ij2idx(i,j))
+    end do
+  end do  
+  
   ! Clean up
-  call tmp_incr%final()
-
+  call field%final()
+  call fieldset%final()
 end subroutine soca_2d_convol
 
 
@@ -398,20 +413,32 @@ end subroutine soca_2d_convol
 !! used by soca_cov::sqrt_C_mult()
 !! \relates soca_covariance_mod::soca_cov
 subroutine soca_2d_sqrt_convol(dx, horiz_convol, geom, pert_scale)
-  real(kind=kind_real), intent(inout) :: dx(:,:)
+  real(kind=kind_real), allocatable, intent(inout) :: dx(:,:)
   type(bump_type),      intent(inout) :: horiz_convol
   type(soca_geom),         intent(in) :: geom
   real(kind=kind_real),    intent(in) :: pert_scale
 
-  type(fieldset_type) :: tmp_incr
+  type(fieldset_type) :: fieldset
+  type(atlas_field)   :: field
+  real(kind_real), pointer :: real_ptr(:,:)
+  integer :: i, j, n
+
   type(atlas_field) :: acv
   integer, parameter :: rseed = 1 ! constant for reproducability of tests
                                   ! TODO: pass seed through config
   integer :: nn
   real(kind=kind_real), pointer :: ptr(:)
 
-  ! Allocate ATLAS tmp_increment and make copy of dx
-  call geom%struct2atlas(dx(:,:), tmp_incr)
+  ! array to atlas fieldset
+  fieldset = atlas_fieldset()
+  field = geom%functionspace%create_field('var', kind=atlas_real(kind_real), levels=1)
+  call fieldset%add(field)
+  call field%data(real_ptr)
+  do j=geom%jsc,geom%jec
+    do i=geom%isc,geom%iec
+      real_ptr(1,geom%atlas_ij2idx(i,j)) = dx(i,j)
+    end do
+  end do  
 
   ! Get control variable size
   call horiz_convol%get_cv_size(nn)
@@ -422,14 +449,19 @@ subroutine soca_2d_sqrt_convol(dx, horiz_convol, geom, pert_scale)
   ptr = pert_scale * ptr
 
   ! Apply C^1/2
-  call horiz_convol%apply_nicas_sqrt(acv, tmp_incr, 0)
+  call horiz_convol%apply_nicas_sqrt(acv, fieldset, 0)
 
-  ! Back to structured grid
-  call geom%atlas2struct(dx(:,:), tmp_incr)
+  ! atlas to array
+  do j=geom%jsc,geom%jec
+    do i=geom%isc,geom%iec
+      dx(i,j) = real_ptr(1,geom%atlas_ij2idx(i,j))
+    end do
+  end do  
 
   ! Clean up
   call acv%final()
-  call tmp_incr%final()
+  call field%final()
+  call fieldset%final()
 
 end subroutine soca_2d_sqrt_convol
 
diff --git a/src/soca/Fields/soca_fields_mod.F90 b/src/soca/Fields/soca_fields_mod.F90
index d531cd9e7..b15030e3b 100644
--- a/src/soca/Fields/soca_fields_mod.F90
+++ b/src/soca/Fields/soca_fields_mod.F90
@@ -1,4 +1,4 @@
-! (C) Copyright 2017-2021 UCAR
+! (C) Copyright 2017-2023 UCAR
 !
 ! This software is licensed under the terms of the Apache Licence Version 2.0
 ! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
@@ -30,7 +30,7 @@ module soca_fields_mod
 use fms_mod,    only: write_data, set_domain
 use MOM_remapping, only : remapping_CS, initialize_remapping, remapping_core_h, &
                           end_remapping
-use mpp_domains_mod, only : mpp_update_domains, mpp_update_domains_ad
+use mpp_domains_mod, only : mpp_update_domains
 
 ! SOCA modules
 use soca_fields_metadata_mod, only : soca_field_metadata
@@ -228,9 +228,6 @@ module soca_fields_mod
   !> copybrief soca_fields_to_fieldset \see soca_fields_to_fieldset
   procedure :: to_fieldset  => soca_fields_to_fieldset
 
-  !> copybrief soca_fields_to_fieldset_ad \see soca_fields_to_fieldset_ad
-  procedure :: to_fieldset_ad  => soca_fields_to_fieldset_ad
-
   procedure :: from_fieldset => soca_fields_from_fieldset
 
   !> \}
@@ -1557,9 +1554,8 @@ subroutine soca_fields_to_fieldset(self, vars, afieldset)
   type(atlas_fieldset), intent(inout) :: afieldset
 
   type(atlas_field) :: afield
-  integer :: v, z
+  integer :: v, n, i, j
   type(soca_field), pointer :: field
-  real(kind=kind_real), pointer :: mask(:,:) => null() !< field mask
   type(atlas_metadata) :: meta
   real(kind=kind_real), pointer :: real_ptr(:,:)
 
@@ -1576,7 +1572,7 @@ subroutine soca_fields_to_fieldset(self, vars, afieldset)
     if (afieldset%has_field(vars%variable(v))) then
       afield = afieldset%field(vars%variable(v))
     else
-      afield = self%geom%functionspaceInchalo%create_field( &
+      afield = self%geom%functionspace%create_field( &
         name=vars%variable(v), kind=atlas_real(kind_real), levels=field%nz)
       meta = afield%metadata()
       call meta%set('interp_type', 'default')
@@ -1588,53 +1584,16 @@ subroutine soca_fields_to_fieldset(self, vars, afieldset)
 
     ! create and fill field
     call afield%data(real_ptr)
-    do z=1,field%nz
-        real_ptr(z,:) = pack(field%val(:,:, z), mask=self%geom%valid_halo_mask)
+    do j=self%geom%jsc,self%geom%jec
+      do i=self%geom%isc,self%geom%iec
+        real_ptr(:, self%geom%atlas_ij2idx(i,j)) = field%val(i,j,:)
+      end do
     end do
-    call afield%final()
-
-  end do
-end subroutine
-
-! ------------------------------------------------------------------------------
-!> Adjoint of get fields used by the interpolation.
-!!
-!! The fields that are input ahave  have halos (minus the invalid and duplicate halo points)
-subroutine soca_fields_to_fieldset_ad(self, vars, afieldset)
-  class(soca_fields),   intent(in) :: self
-  type(oops_variables), intent(in) :: vars
-  type(atlas_fieldset), intent(in) :: afieldset
-
-  integer :: v, z
-  integer :: is, ie, js, je
-  type(soca_field), pointer :: field
-  type(atlas_field) :: afield
-  real(kind=kind_real), pointer :: real_ptr(:,:)
-  real(kind=kind_real), pointer :: tmp(:,:)
-
-  ! start/stop idx, assuming halo
-  is = self%geom%isd; ie = self%geom%ied
-  js = self%geom%jsd; je = self%geom%jed
 
-  allocate(tmp(is:ie, js:je))
-
-  do v=1,vars%nvars()
-    call self%get(vars%variable(v), field)
-    afield = afieldset%field(vars%variable(v))
-
-    tmp = 0.0
-    call afield%data(real_ptr)
-    do z=1,field%nz
-      tmp = unpack(real_ptr(z,:), self%geom%valid_halo_mask, 0.0_kind_real)
-      call mpp_update_domains_ad(tmp, self%geom%Domain%mpp_domain, complete=.true.)
-      field%val(:,:,z) = field%val(:,:,z) + tmp
-    end do
     call afield%final()
-
   end do
 end subroutine
 
-
 ! ------------------------------------------------------------------------------
 !> Set the our values from an atlas fieldset
 subroutine soca_fields_from_fieldset(self, vars, afieldset)
@@ -1642,7 +1601,7 @@ subroutine soca_fields_from_fieldset(self, vars, afieldset)
   type(oops_variables),       intent(in)    :: vars
   type(atlas_fieldset),       intent(in)    :: afieldset
 
-  integer :: jvar, i, jz
+  integer :: jvar, i, j, n, f
   real(kind=kind_real), pointer :: real_ptr(:,:)
   logical :: var_found
   character(len=1024) :: fieldname
@@ -1654,20 +1613,20 @@ subroutine soca_fields_from_fieldset(self, vars, afieldset)
 
   do jvar = 1,vars%nvars()
     var_found = .false.
-    do i=1,size(self%fields)
-      field => self%fields(i)
+    do f=1,size(self%fields)
+      field => self%fields(f)
       if (trim(vars%variable(jvar))==trim(field%name)) then
         ! Get field
         afield = afieldset%field(vars%variable(jvar))
 
         ! Copy data
         call afield%data(real_ptr)
-        do jz=1,field%nz
-          ! NOTE, any missing values from unpacking are filled with 0.0
-          field%val(:,:,jz) = unpack(real_ptr(jz,:), self%geom%valid_halo_mask, 0.0_kind_real)
+        do j=self%geom%jsc,self%geom%jec
+          do i=self%geom%isc,self%geom%iec
+            field%val(i,j,:) = real_ptr(:, self%geom%atlas_ij2idx(i,j))
+          end do
         end do
-
-        ! Release pointer
+        
         call afield%final()
 
         ! Set flag
@@ -1677,7 +1636,6 @@ subroutine soca_fields_from_fieldset(self, vars, afieldset)
     end do
     if (.not.var_found) call abor1_ftn('variable '//trim(vars%variable(jvar))//' not found in increment')
   end do
-
 end subroutine
 
 end module soca_fields_mod
diff --git a/src/soca/Geometry/Geometry.cc b/src/soca/Geometry/Geometry.cc
index cc344038b..2873837c0 100644
--- a/src/soca/Geometry/Geometry.cc
+++ b/src/soca/Geometry/Geometry.cc
@@ -1,17 +1,23 @@
 /*
- * (C) Copyright 2017-2021 UCAR
+ * (C) Copyright 2017-2023 UCAR
  *
  * This software is licensed under the terms of the Apache Licence Version 2.0
  * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
  */
 
-#include "atlas/grid.h"
-#include "atlas/util/Config.h"
+#include <algorithm>
 
-#include "eckit/config/YAMLConfiguration.h"
+#include "atlas/functionspace.h"
+#include "atlas/mesh/actions/BuildHalo.h"
+#include "atlas/mesh/Mesh.h"
+#include "atlas/mesh/MeshBuilder.h"
+#include "atlas/output/Gmsh.h"
+
+#include "eckit/config/Configuration.h"
 
 #include "soca/Geometry/Geometry.h"
 
+
 // -----------------------------------------------------------------------------
 namespace soca {
 
@@ -27,27 +33,90 @@ namespace soca {
 
     soca_geo_setup_f90(keyGeom_, &conf, &comm);
 
-    // Set ATLAS lonlat and function space (with and without halos)
-    atlas::FieldSet lonlat;
-    soca_geo_lonlat_f90(keyGeom_, lonlat.get());
-    functionSpace_ = atlas::functionspace::PointCloud(lonlat->field("lonlat"));
-    functionSpaceIncHalo_ = atlas::functionspace::PointCloud(lonlat->field("lonlat_inc_halos"));
-
-    // Set ATLAS function space pointer in Fortran
-    soca_geo_set_atlas_functionspace_pointer_f90(
-      keyGeom_, functionSpace_.get(), functionSpaceIncHalo_.get());
-
-    // Fill ATLAS fieldset
-    soca_geo_to_fieldset_f90(keyGeom_, fields_.get());
-
-    // messy, fix this
     // generate the grid ONLY if being run under the gridgen application.
-    // also, if true, then don't bother with the kdtree generation in the next step.
     if (gen) {
       soca_geo_gridgen_f90(keyGeom_);
-      return;
     }
+
+    // setup the atlas functionspace
+    {
+      using atlas::gidx_t;
+      using atlas::idx_t;
+
+      // get the number of nodes and cells owned by this PE
+      int num_nodes;
+      int num_quad_elements;
+      soca_geo_get_mesh_size_f90(keyGeom_, num_nodes, num_quad_elements);
+
+      // get the mesh connectivity from the soca fortran
+      std::vector<double> lons(num_nodes);
+      std::vector<double> lats(num_nodes);
+      std::vector<int> ghosts(num_nodes);
+      std::vector<int> global_indices(num_nodes);
+      std::vector<int> remote_indices(num_nodes);
+      std::vector<int> partitions(num_nodes);
+      const int num_quad_nodes = num_quad_elements * 4;
+      std::vector<int> raw_quad_nodes(num_quad_nodes);
+      soca_geo_gen_mesh_f90(keyGeom_,
+        num_nodes, lons.data(), lats.data(), ghosts.data(), global_indices.data(),
+        remote_indices.data(), partitions.data(),
+        num_quad_nodes, raw_quad_nodes.data());
+
+      // calculate per-PE global quad numbering offset
+      std::vector<int> num_elements_per_rank(comm_.size());
+      comm_.allGather(num_quad_elements, num_elements_per_rank.begin(),
+                      num_elements_per_rank.end());
+      int global_element_index = 0;
+      for (size_t i = 0; i < comm_.rank(); ++i) {
+        global_element_index += num_elements_per_rank[i];
+      }
+
+      // convert some of the temporary arrays into a form atlas expects
+      std::vector<gidx_t> atlas_global_indices(num_nodes);
+      std::transform(global_indices.begin(), global_indices.end(), atlas_global_indices.begin(),
+        [](const int index) {return atlas::gidx_t{index};});
+      std::vector<idx_t> atlas_remote_indices(num_nodes);
+      std::transform(remote_indices.begin(), remote_indices.end(), atlas_remote_indices.begin(),
+        [](const int index) {return atlas::idx_t{index};});
+      std::vector<std::array<gidx_t, 3>> tri_boundary_nodes{};  // MOM does not have triangles
+      std::vector<gidx_t> tri_global_indices{};  // MOM does not have triangles
+      std::vector<std::array<gidx_t, 4>> quad_boundary_nodes(num_quad_elements);
+      std::vector<gidx_t> quad_global_indices(num_quad_elements);
+      for (size_t quad = 0; quad < num_quad_elements; ++quad) {
+        for (size_t i = 0; i < 4; ++i) {
+          quad_boundary_nodes[quad][i] = raw_quad_nodes[4*quad + i];
+        }
+        quad_global_indices[quad] = global_element_index++;
+      }
+
+      // build the mesh!
+      const atlas::idx_t remote_index_base = 1;  // 1-based indexing from Fortran
+      eckit::LocalConfiguration config{};
+      config.set("mpi_comm", comm_.name());
+      const atlas::mesh::MeshBuilder mesh_builder{};
+      atlas::Mesh mesh = mesh_builder(
+        lons, lats, ghosts,
+        atlas_global_indices, atlas_remote_indices, remote_index_base, partitions,
+        tri_boundary_nodes, tri_global_indices,
+        quad_boundary_nodes, quad_global_indices, config);
+      atlas::mesh::actions::build_halo(mesh, 1);
+      functionSpace_ = atlas::functionspace::NodeColumns(mesh, config);
+
+      // optionally save output for viewing with gmsh
+      if (conf.getBool("gmsh save", false)) {
+        std::string filename = conf.getString("gmsh filename", "out.msh");
+        atlas::output::Gmsh gmsh(filename,
+            atlas::util::Config("coordinates", "xyz")
+            | atlas::util::Config("ghost", true));  // enables viewing halos per task
+        gmsh.write(mesh);
+      }
+    }
+
+    // Set ATLAS function space in Fortran, and fill in the
+    // geometry fieldset from the fortran side.
+    soca_geo_init_atlas_f90(keyGeom_, functionSpace_.get(), fields_.get());
   }
+
   // -----------------------------------------------------------------------------
   Geometry::Geometry(const Geometry & other)
     : comm_(other.comm_),
@@ -56,16 +125,8 @@ namespace soca {
     const int key_geo = other.keyGeom_;
     soca_geo_clone_f90(keyGeom_, key_geo);
 
-    functionSpace_ = atlas::functionspace::PointCloud(other.functionSpace_->lonlat());
-    functionSpaceIncHalo_ = atlas::functionspace::PointCloud(other.functionSpaceIncHalo_->lonlat());
-    soca_geo_set_atlas_functionspace_pointer_f90(keyGeom_, functionSpace_.get(),
-                                                 functionSpaceIncHalo_.get());
-
-    fields_ = atlas::FieldSet();
-    for (int jfield = 0; jfield < other.fields_->size(); ++jfield) {
-      atlas::Field atlasField = other.fields_->field(jfield);
-      fields_->add(atlasField);
-    }
+    functionSpace_ = atlas::functionspace::NodeColumns(other.functionSpace_);
+    soca_geo_init_atlas_f90(keyGeom_, functionSpace_.get(), fields_.get());
   }
   // -----------------------------------------------------------------------------
   Geometry::~Geometry() {
@@ -110,14 +171,32 @@ namespace soca {
   // -----------------------------------------------------------------------------
   void Geometry::latlon(std::vector<double> & lats, std::vector<double> & lons,
       const bool halo) const {
-    // get the number of gridpoints
-    int gridSize;
-    soca_geo_gridsize_f90(keyGeom_, halo, gridSize);
+    // get the number of total grid points (including halo)
+    int gridSizeWithHalo = functionSpace_.size();
+    auto vLonlat = atlas::array::make_view<double, 2>(functionSpace_.lonlat());
+
+    // count the number of owned non-ghost points (isn't there an atlas function for this??)
+    auto vGhost = atlas::array::make_view<int, 1>(functionSpace_.ghost());
+    int gridSizeNoHalo = 0;
+    for (size_t i = 0; i < gridSizeWithHalo; i++) {
+      if (vGhost(i) == 0) gridSizeNoHalo++;
+    }
 
-    // get the lat/lon of those gridpoints
-    lats.resize(gridSize);
+    // allocate arrays
+    int gridSize = (halo) ? gridSizeWithHalo : gridSizeNoHalo;
     lons.resize(gridSize);
-    soca_geo_gridlatlon_f90(keyGeom_, halo, gridSize, lats.data(), lons.data());
+    lats.resize(gridSize);
+
+    // fill
+    int idx = 0;
+    for (size_t i=0; i < gridSizeWithHalo; i++) {
+      if (!halo && vGhost(i)) continue;
+      double lon = vLonlat(i, 0);
+      double lat = vLonlat(i, 1);
+      lats[idx] = lat;
+      lons[idx++] = lon;
+    }
+    ASSERT(idx == gridSize);
   }
 
 }  // namespace soca
diff --git a/src/soca/Geometry/Geometry.h b/src/soca/Geometry/Geometry.h
index 2f4ab989b..17bcb1320 100644
--- a/src/soca/Geometry/Geometry.h
+++ b/src/soca/Geometry/Geometry.h
@@ -69,8 +69,8 @@ namespace soca {
       const int& toFortran() const {return keyGeom_;}
       const eckit::mpi::Comm & getComm() const {return comm_;}
 
-      const atlas::FunctionSpace & functionSpace() const {return functionSpaceIncHalo_;}
-      atlas::FunctionSpace & functionSpace() {return functionSpaceIncHalo_;}
+      const atlas::FunctionSpace & functionSpace() const {return functionSpace_;}
+      atlas::FunctionSpace & functionSpace() {return functionSpace_;}
       const atlas::FieldSet & fields() const {return fields_;}
       atlas::FieldSet & fields() {return fields_;}
 
@@ -84,8 +84,7 @@ namespace soca {
       int keyGeom_;
       const eckit::mpi::Comm & comm_;
       FmsInput fmsinput_;
-      atlas::FunctionSpace functionSpace_;
-      atlas::FunctionSpace functionSpaceIncHalo_;
+      atlas::functionspace::NodeColumns functionSpace_;
       atlas::FieldSet fields_;
   };
   // -----------------------------------------------------------------------------
diff --git a/src/soca/Geometry/GeometryFortran.h b/src/soca/Geometry/GeometryFortran.h
index 77b465a81..1d4d855f5 100644
--- a/src/soca/Geometry/GeometryFortran.h
+++ b/src/soca/Geometry/GeometryFortran.h
@@ -16,6 +16,7 @@
 namespace atlas {
   namespace field {
     class FieldSetImpl;
+    class FieldImpl;
   }
   namespace functionspace {
     class FunctionSpaceImpl;
@@ -32,14 +33,14 @@ namespace soca {
                             const eckit::Configuration * const &,
                             const eckit::mpi::Comm *);
 
-    void soca_geo_lonlat_f90(const F90geom &,
-                             atlas::field::FieldSetImpl *);
-    void soca_geo_set_atlas_functionspace_pointer_f90(const F90geom &,
+    void soca_geo_init_atlas_f90(const F90geom &,
                       atlas::functionspace::FunctionSpaceImpl *,
-                      atlas::functionspace::FunctionSpaceImpl *);
-    void soca_geo_to_fieldset_f90(const F90geom &,
-                                  atlas::field::FieldSetImpl *);
-
+                      atlas::field::FieldSetImpl *);
+    void soca_geo_get_mesh_size_f90(const F90geom &, int &, int&);
+    void soca_geo_gen_mesh_f90(
+      const F90geom &,
+      const int &, double[], double[], int[], int[], int[], int[],
+      const int &, int[]);
     void soca_geo_clone_f90(F90geom &, const F90geom &);
     void soca_geo_gridgen_f90(const F90geom &);
     void soca_geo_delete_f90(F90geom &);
@@ -48,9 +49,6 @@ namespace soca {
     void soca_geo_get_num_levels_f90(const F90geom &, const oops::Variables &,
                                     const size_t &, size_t[]);
     void soca_geo_iterator_dimension_f90(const F90geom &, int &);
-
-    void soca_geo_gridsize_f90(const F90geom &, const bool &, int &);
-    void soca_geo_gridlatlon_f90(const F90geom &, const bool &, const int &, double[], double[]);
   }
 }  // namespace soca
 #endif  // SOCA_GEOMETRY_GEOMETRYFORTRAN_H_
diff --git a/src/soca/Geometry/soca_geom.interface.F90 b/src/soca/Geometry/soca_geom.interface.F90
index 7f440037d..6157a413c 100644
--- a/src/soca/Geometry/soca_geom.interface.F90
+++ b/src/soca/Geometry/soca_geom.interface.F90
@@ -1,4 +1,4 @@
-! (C) Copyright 2017-2021 UCAR
+! (C) Copyright 2017-2023UCAR
 !
 ! This software is licensed under the terms of the Apache Licence Version 2.0
 ! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
@@ -7,13 +7,15 @@
 !> C++ interfaces for soca_geom_mod::soca_geom
 module soca_geom_mod_c
 
-use atlas_module, only: atlas_fieldset, atlas_functionspace_pointcloud
+use atlas_module, only: atlas_fieldset, atlas_functionspace_NodeColumns, atlas_field
 use fckit_configuration_module, only: fckit_configuration
 use fckit_mpi_module,           only: fckit_mpi_comm
 use iso_c_binding
 use oops_variables_mod, only: oops_variables
 use kinds, only: kind_real
 
+use mpp_domains_mod, only : mpp_update_domains
+
 ! soca modules
 use soca_geom_mod, only: soca_geom
 use soca_fields_metadata_mod, only : soca_field_metadata
@@ -57,53 +59,23 @@ end subroutine soca_geo_setup_c
 
 
 ! --------------------------------------------------------------------------------------------------
-!> C++ interface for soca_geom_mod::soca_geom::set_atlas_lonlat()
-subroutine soca_geo_lonlat_c(c_key_self, c_afieldset)  bind(c,name='soca_geo_lonlat_f90')
-  integer(c_int), intent(in) :: c_key_self
-  type(c_ptr), intent(in), value :: c_afieldset
-
-  type(soca_geom), pointer :: self
-  type(atlas_fieldset) :: afieldset
-
-  call soca_geom_registry%get(c_key_self,self)
-  afieldset = atlas_fieldset(c_afieldset)
-
-  call self%lonlat(afieldset)
-end subroutine soca_geo_lonlat_c
-
-
-! --------------------------------------------------------------------------------------------------
-!> C++ interface to get atlas functionspace pointr from  soca_geom_mod::soca_geom
-subroutine soca_geo_set_atlas_functionspace_pointer_c(c_key_self, c_functionspace, c_functionspaceIncHalo) &
-  bind(c,name='soca_geo_set_atlas_functionspace_pointer_f90')
+!> C++ interface to get atlas functionspace pointr from  soca_geom_mod::soca_geom,
+! and then fill in the atlas fieldset with the required geometry fields
+subroutine soca_geo_init_atlas_c(c_key_self, c_functionspace, c_fieldset) &
+    bind(c,name='soca_geo_init_atlas_f90')
   integer(c_int), intent(in)     :: c_key_self
-  type(c_ptr), intent(in), value :: c_functionspace, c_functionspaceIncHalo
+  type(c_ptr), intent(in), value :: c_functionspace, c_fieldset
 
   type(soca_geom),pointer :: self
 
   call soca_geom_registry%get(c_key_self,self)
+  self%functionspace = atlas_functionspace_NodeColumns(c_functionspace)
 
-  self%functionspace = atlas_functionspace_pointcloud(c_functionspace)
-  self%functionspaceIncHalo = atlas_functionspace_pointcloud(c_functionspaceIncHalo)
-end subroutine soca_geo_set_atlas_functionspace_pointer_c
-
-
-! --------------------------------------------------------------------------------------------------
-!> C++ interface for soca_geom_mod::soca_geom::fill_atlas_fieldset()
-subroutine soca_geo_to_fieldset_c(c_key_self, c_afieldset) &
- & bind(c,name='soca_geo_to_fieldset_f90')
-
-  integer(c_int),     intent(in) :: c_key_self
-  type(c_ptr), value, intent(in) :: c_afieldset
-
-  type(soca_geom), pointer :: self
-  type(atlas_fieldset) :: afieldset
-
-  call soca_geom_registry%get(c_key_self,self)
-  afieldset = atlas_fieldset(c_afieldset)
-
-  call self%to_fieldset(afieldset)
-end subroutine soca_geo_to_fieldset_c
+  ! fill in the geometry fieldset  
+  self%fieldset = atlas_fieldset(c_fieldset)
+  call self%init_fieldset()
+  
+end subroutine soca_geo_init_atlas_c
 
 
 ! ------------------------------------------------------------------------------
@@ -215,56 +187,115 @@ subroutine soca_geo_get_num_levels_c(c_key_self, c_vars, c_levels_size, c_levels
 end subroutine soca_geo_get_num_levels_c
 
 ! ------------------------------------------------------------------------------
-!> Get the number of points valid in the local grid (skipping masked points)
-subroutine soca_geo_gridsize_c(c_key_self, c_halo, c_size) &
-           bind(c, name='soca_geo_gridsize_f90')
-  integer(c_int),    intent(in) :: c_key_self
-  logical(c_bool),   intent(in) :: c_halo    !< true if halo should be included in number of points
-  integer(c_int),    intent(out):: c_size    !< the resulting number of gridpoints
+! Get the number of nodes and number of quadrilaterals that 
+! will be generated by a subsequent call to soca_geo_gen_mesh
+subroutine soca_geo_get_mesh_size_c(c_key_self, c_nodes, c_quads) bind(c, name='soca_geo_get_mesh_size_f90')
+  integer(c_int), intent(in)  :: c_key_self
+  integer(c_int), intent(out) :: c_nodes, c_quads
 
   type(soca_geom), pointer :: self
+  logical, allocatable :: valid_nodes(:,:), valid_cells(:,:)
 
   call soca_geom_registry%get(c_key_self, self)
-  if (c_halo) then
-    c_size = self%ngrid_halo_valid
-  else
-    c_size = self%ngrid
-  end if
-end subroutine soca_geo_gridsize_c
-
-! ------------------------------------------------------------------------------
-!> Get the lat/lons for the h grid
-!> If c_halo is true, then only the "valid" halo points are used
-!> (i.e. duplicate halo points and invalid halo points are skipped)
-subroutine soca_geo_gridlatlon_c(c_key_self, c_halo, c_size, &
-    c_lat, c_lon) bind(c, name='soca_geo_gridlatlon_f90')
+  
+  call self%mesh_valid_nodes_cells(valid_nodes, valid_cells)
+  c_nodes = count(valid_nodes)
+  c_quads = count(valid_cells)
 
-  integer(c_int),    intent(in) :: c_key_self
-  logical(c_bool),   intent(in) :: c_halo
-  integer(c_int),    intent(in) :: c_size
-  real(c_double),    intent(inout)  :: c_lat(c_size), c_lon(c_size)
+end subroutine
 
+! ------------------------------------------------------------------------------
+! Generate the node and quadrilateral information that is needed
+! by the C++ side of Geometry::Geometry() to create a connected mesh
+! in atlas.
+subroutine soca_geo_gen_mesh_c(c_key_self, c_nodes, c_lon, c_lat, c_ghosts, c_global_idx, c_remote_idx, c_partition, &
+    c_quad_nodes, c_quad_node_list) bind(c, name='soca_geo_gen_mesh_f90')
+  integer(c_int), intent(in)    :: c_key_self
+  integer(c_int), intent(in)    :: c_nodes, c_quad_nodes
+  integer(c_int), intent(inout) :: c_ghosts(c_nodes), c_global_idx(c_nodes), c_remote_idx(c_nodes), c_partition(c_nodes)
+  integer(c_int), intent(inout) :: c_quad_node_list(c_quad_nodes)
+  real(c_double), intent(inout) :: c_lon(c_nodes), c_lat(c_nodes)
+
+  integer :: idx, i, j
+  integer :: nx, ny
+  integer, allocatable :: global_idx(:,:), local_idx(:,:), partition(:,:)
+  logical, allocatable :: valid_nodes(:,:), valid_cells(:,:)
+  
   type(soca_geom), pointer :: self
-  real(kind=kind_real),     pointer :: lon(:,:) => null()  !< field lon
-  real(kind=kind_real),     pointer :: lat(:,:) => null()  !< field lat
 
   call soca_geom_registry%get(c_key_self, self)
-  
-  ! TODO: re-enable the ability to get lat/lon for different grids?
-  ! For now just use h grid
-  lon => self%lon
-  lat => self%lat
-  
-  if (c_halo) then
-    ! get all lat/lon points, except for duplicate or invalid halo points
-    c_lat(:) = pack(self%lat, mask=self%valid_halo_mask)
-    c_lon(:) = pack(self%lon, mask=self%valid_halo_mask)
-  else
-    ! get all non-halo lat/lon points
-    c_lat(:) = pack(self%lat(self%isc:self%iec, self%jsc:self%jec), .true.)
-    c_lon(:) = pack(self%lon(self%isc:self%iec, self%jsc:self%jec), .true.)
+
+  ! parameters pulled from grid
+  nx=self%domain%NIGLOBAL
+  ny=self%domain%NJGLOBAL
+
+  ! find global/local indexes, and the PE owner of each grid point
+  ! (requiring a halo communications)
+  allocate(global_idx(self%isd:self%ied, self%jsd:self%jed))
+  allocate(local_idx(self%isd:self%ied, self%jsd:self%jed))
+  allocate(partition(self%isd:self%ied, self%jsd:self%jed))
+  idx=0 ! local index
+  do j=self%jsc,self%jec
+    do i=self%isc, self%iec
+      idx = idx + 1
+      local_idx(i,j) = idx
+      global_idx(i,j) = (j-1)*nx + i
+      partition(i,j) = self%f_comm%rank()
+    end do
+  end do
+  call mpp_update_domains(local_idx, self%Domain%mpp_domain)
+  call mpp_update_domains(global_idx, self%Domain%mpp_domain)
+  call mpp_update_domains(partition, self%Domain%mpp_domain)
+
+  ! find which quads / vertices are we going to skip (in case of non cyclic or tripolar_fold special cases)
+  call self%mesh_valid_nodes_cells(valid_nodes, valid_cells)
+
+
+  ! fill in the node arrays to return to C++.
+  ! AND while were at it, save the atlas array to 
+  ! soca fortran array mapping for use later by to_fieldset
+  ! and from_fieldset
+  c_ghosts = 1
+  idx=1
+  do j=self%jsc,self%jec+1
+    do i=self%isc, self%iec+1
+      if (.not. valid_nodes(i,j)) cycle
+
+      ! fill in the node information arrays
+      c_lon(idx) = self%lon(i,j)
+      c_lat(idx) = self%lat(i,j)
+      if(j .le. self%jec .and. i .le. self%iec) then
+        c_ghosts(idx) = 0
+        self%atlas_ij2idx(i,j) = idx
+      end if
+      c_global_idx(idx) = global_idx(i,j)
+      c_remote_idx(idx) = local_idx(i,j)
+      c_partition(idx) = partition(i,j)
+
+      idx = idx + 1
+    end do
+  end do
+  if (c_nodes /= idx-1) then
+    call abor1_ftn("ERROR: c_nodes != idx-1")
   end if
-end subroutine
 
+  ! fill in the quad node list arrays
+  idx=1
+  do j=self%jsc,self%jec
+    do i=self%isc, self%iec
+      if(.not. valid_cells(i,j)) cycle
+
+      c_quad_node_list(idx  ) = global_idx(i  ,j  )
+      c_quad_node_list(idx+1) = global_idx(i  ,j+1)
+      c_quad_node_list(idx+2) = global_idx(i+1,j+1)
+      c_quad_node_list(idx+3) = global_idx(i+1,j  )
+      
+      idx = idx + 4
+    end do
+  end do
+  if (c_quad_nodes /= idx-1) then
+    call abor1_ftn("ERROR: c_quad_nodes != idx-1")
+  end if
+end subroutine
 ! ------------------------------------------------------------------------------
 end module soca_geom_mod_c
diff --git a/src/soca/Geometry/soca_geom_mod.F90 b/src/soca/Geometry/soca_geom_mod.F90
index 52c1e1eec..aa19445a0 100644
--- a/src/soca/Geometry/soca_geom_mod.F90
+++ b/src/soca/Geometry/soca_geom_mod.F90
@@ -1,4 +1,4 @@
-! (C) Copyright 2017-2021 UCAR
+! (C) Copyright 2017-2023 UCAR
 !
 ! This software is licensed under the terms of the Apache Licence Version 2.0
 ! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
@@ -7,8 +7,8 @@
 !> Geometry module
 module soca_geom_mod
 
-! jedi modules
-use atlas_module, only: atlas_functionspace_pointcloud, atlas_fieldset, &
+  ! jedi modules
+use atlas_module, only: atlas_functionspace_NodeColumns, atlas_fieldset, &
     atlas_field, atlas_real, atlas_integer, atlas_geometry, atlas_indexkdtree
 use fckit_configuration_module, only: fckit_configuration
 use fckit_mpi_module, only: fckit_mpi_comm
@@ -24,7 +24,8 @@ module soca_geom_mod
 use MOM_domains, only : MOM_domain_type
 use MOM_EOS,         only : EOS_type
 use mpp_domains_mod, only : mpp_get_compute_domain, mpp_get_data_domain, &
-                            mpp_get_global_domain, mpp_update_domains
+                            mpp_get_global_domain, mpp_update_domains, &
+                            CYCLIC_GLOBAL_DOMAIN, FOLD_NORTH_EDGE                            
 ! soca modules
 use soca_fields_metadata_mod, only : soca_fields_metadata
 use soca_mom6, only: soca_mom6_config, soca_mom6_init, soca_geomdomain_init
@@ -80,11 +81,6 @@ module soca_geom_mod
     integer :: iterator_dimension
     !> \}
 
-    integer :: ngrid !< number of non-halo gridpoints
-    integer :: ngrid_halo !< number of gridpoints, including halo
-    integer :: ngrid_halo_valid !< number of gridpoints, including halo, but
-                                !! excluding duplicate or invalid halo points
-
     !> \name grid latitude/longitude
     !! \{
     real(kind=kind_real), allocatable, dimension(:)   :: lonh !< cell center nominal longitude
@@ -121,7 +117,6 @@ module soca_geom_mod
     real(kind=kind_real), allocatable, dimension(:,:) :: distance_from_coast !< distance to closest land grid point (m)
     real(kind=kind_real), allocatable, dimension(:,:,:) :: h !< layer thickness (m)
     real(kind=kind_real), allocatable, dimension(:,:,:) :: h_zstar
-    logical,              allocatable, dimension(:,:) :: valid_halo_mask !< true unless gridpoint is a duplicate or invalid halo
     !> \}
 
     !> instance of the metadata that is read in from a config file upon initialization
@@ -132,9 +127,11 @@ module soca_geom_mod
     character(len=:), allocatable :: rossby_file !< filename of rossby radius input file (if used)
 
     type(fckit_mpi_comm) :: f_comm !< MPI communicator
-    type(atlas_functionspace_pointcloud) :: functionspace
-    type(atlas_functionspace_pointcloud) :: functionspaceInchalo
 
+    !> mesh parameters
+    type(atlas_functionspace_NodeColumns) :: functionspace
+    integer, allocatable :: atlas_ij2idx(:,:)
+    type(atlas_fieldset) :: fieldset !< the geom fields (area, mask, etc)
 
   contains
 
@@ -145,11 +142,8 @@ module soca_geom_mod
     !> \copybrief soca_geom_end \see soca_geom_end
     procedure :: end => soca_geom_end
 
-    !> \copybrief soca_geom_set_atlas_lonlat \see soca_geom_set_atlas_lonlat
-    procedure :: lonlat => soca_geom_lonlat
-
     !> \copybrief soca_geom_fill_atlas_fieldset \see soca_geom_fill_atlas_fieldset
-    procedure :: to_fieldset => soca_geom_to_fieldset
+    procedure :: init_fieldset => soca_geom_init_fieldset
 
     !> \copybrief soca_geom_clone \see soca_geom_clone
     procedure :: clone => soca_geom_clone
@@ -160,15 +154,11 @@ module soca_geom_mod
     !> \copybrief soca_geom_thickness2depth \see soca_geom_thickness2depth
     procedure :: thickness2depth => soca_geom_thickness2depth
 
-    !> \copybrief soca_geom_struct2atlas \see soca_geom_struct2atlas
-    procedure :: struct2atlas => soca_geom_struct2atlas
-
-    !> \copybrief soca_geom_atlas2struct \ see soca_geom_atlas2struct
-    procedure :: atlas2struct => soca_geom_atlas2struct
-
     !> \copybrief soca_geom_write \see soca_geom_write
     procedure :: write => soca_geom_write
 
+    !> \copybrief soca_geom_mesh_valid_nodes_cells \see soca_geom_mesh_valid_nodes_cells
+    procedure :: mesh_valid_nodes_cells => soca_geom_mesh_valid_nodes_cells
 end type soca_geom
 
 ! ------------------------------------------------------------------------------
@@ -227,9 +217,6 @@ subroutine soca_geom_init(self, f_conf, f_comm)
   call mpp_update_domains(self%rossby_radius, self%Domain%mpp_domain)
   call mpp_update_domains(self%distance_from_coast, self%Domain%mpp_domain)
 
-  ! calculate which halo points are duplicates / invalid
-  call soca_geom_find_invalid_halo(self)
-
   ! Set output option for local geometry
   if ( .not. f_conf%get("save_local_domain", self%save_local_domain) ) &
      self%save_local_domain = .false.
@@ -273,102 +260,81 @@ subroutine soca_geom_end(self)
   if (allocated(self%distance_from_coast)) deallocate(self%distance_from_coast)
   if (allocated(self%h))             deallocate(self%h)
   if (allocated(self%h_zstar))       deallocate(self%h_zstar)
-  if (allocated(self%valid_halo_mask)) deallocate(self%valid_halo_mask)
   nullify(self%Domain)
   call self%functionspace%final()
-  call self%functionspaceIncHalo%final()
 
 end subroutine soca_geom_end
 
 
 ! --------------------------------------------------------------------------------------------------
-!> Set ATLAS lonlat fieldset
+!> Fill ATLAS fieldset
 !!
 !! \related soca_geom_mod::soca_geom
-subroutine soca_geom_lonlat(self, afieldset)
+subroutine soca_geom_init_fieldset(self)
   class(soca_geom),  intent(inout) :: self
-  type(atlas_fieldset), intent(inout) :: afieldset
-
-  real(kind_real), pointer :: real_ptr(:,:)
-  type(atlas_field) :: afield
 
-  ! Create lon/lat field
-  afield = atlas_field(name="lonlat", kind=atlas_real(kind_real), shape=(/2,self%ngrid/))
-  call afield%data(real_ptr)
-  real_ptr(1,:) = pack(self%lon(self%isc:self%iec,self%jsc:self%jec),.true.)
-  real_ptr(2,:) = pack(self%lat(self%isc:self%iec,self%jsc:self%jec),.true.)
-  call afieldset%add(afield)
-  call afield%final()
+  integer :: i, j, n, jz
+  type(atlas_field) :: fArea, fInterpMask, fVertCoord, fGmask, fOwned, fRossby
+  real(kind=kind_real), pointer :: vArea(:,:), vInterpMask(:,:), vVertCoord(:,:), vRossby(:,:)
+  integer, pointer :: vGmask(:,:), vOwned(:,:)
 
-  afield = atlas_field(name="lonlat_inc_halos", kind=atlas_real(kind_real), shape=(/2,self%ngrid_halo_valid/))
-  call afield%data(real_ptr)
-  real_ptr(1,:) = pack(self%lon, self%valid_halo_mask)
-  real_ptr(2,:) = pack(self%lat, self%valid_halo_mask)
-  call afieldset%add(afield)
-  call afield%final()
-
-end subroutine soca_geom_lonlat
-
-
-! --------------------------------------------------------------------------------------------------
-!> Fill ATLAS fieldset
-!!
-!! \related soca_geom_mod::soca_geom
-subroutine soca_geom_to_fieldset(self, afieldset)
-  class(soca_geom),  intent(inout) :: self
-  type(atlas_fieldset), intent(inout) :: afieldset
-
-  integer :: i, jz, ngrid
-  integer, pointer :: int_ptr(:,:)
-  real(kind=kind_real), pointer :: real_ptr(:,:)
-  integer, allocatable :: hmask(:,:)
-  type(atlas_field) :: afield
-
-  ngrid = (self%ied-self%isd+1)*(self%jed-self%jsd+1)
-
-  ! Add area
-  afield = self%functionspaceInchalo%create_field(name='area', kind=atlas_real(kind_real), levels=1)
-  call afield%data(real_ptr)
-  real_ptr(1,:) = pack(self%cell_area, self%valid_halo_mask)
-  call afieldset%add(afield)
-  call afield%final()
-
-  ! Add vertical unit
-  afield = self%functionspaceInchalo%create_field(name='vert_coord', kind=atlas_real(kind_real), levels=self%nzo)
-  call afield%data(real_ptr)
-  do jz=1,self%nzo
-    real_ptr(jz,:) = real(jz, kind_real)
+  ! create fields, get pointers to their data
+  fArea = self%functionspace%create_field(name='area', kind=atlas_real(kind_real), levels=1)
+  call self%fieldset%add(fArea)
+  call fArea%data(vArea)
+  
+  fInterpMask = self%functionspace%create_field(name='interp_mask', kind=atlas_real(kind_real), levels=1)
+  call self%fieldset%add(fInterpMask)    
+  call fInterpMask%data(vInterpMask)
+
+  fVertCoord = self%functionspace%create_field(name='vert_coord', kind=atlas_real(kind_real), levels=self%nzo)
+  call self%fieldset%add(fVertCoord)
+  call fVertCoord%data(vVertCoord) 
+
+  fGmask = self%functionspace%create_field(name='gmask', kind=atlas_integer(kind(0)), levels=self%nzo)
+  call self%fieldset%add(fGmask)
+  call fGmask%data(vGmask) 
+
+  fOwned = self%functionspace%create_field(name='owned', kind=atlas_integer(kind(0)), levels=1)
+  call self%fieldset%add(fOwned)
+  call fOwned%data(vOwned)
+
+  fRossby = self%functionspace%create_field(name='rossby_radius', kind=atlas_real(kind_real), levels=1)
+  call self%fieldset%add(fRossby)
+  call fRossby%data(vRossby)
+
+  ! set the data
+  vOwned = 0 ! need to set to 0, it's the only one that doesn't get halo update
+  do j=self%jsc,self%jec
+    do i=self%isc,self%iec
+      n = self%atlas_ij2idx(i,j)
+      vArea(1,n) = self%cell_area(i,j)
+      vInterpMask(1,n) = self%mask2d(i,j)
+      vGmask(:, n) = int(self%mask2d(i,j))
+      vOwned(1, n) = 1
+      vRossby(1, n) = self%rossby_radius(i,j)
+    end do
   end do
-  call afieldset%add(afield)
-  call afield%final()
-
-  ! Add geographical mask
-  afield = self%functionspaceInchalo%create_field(name='gmask', kind=atlas_integer(kind(0)), levels=self%nzo)
-  call afield%data(int_ptr)
   do jz=1,self%nzo
-    int_ptr(jz,:) = int(pack(self%mask2d, self%valid_halo_mask))
+    vVertCoord(jz,:) = real(jz, kind_real)
   end do
-  call afieldset%add(afield)
-  call afield%final()
 
-  ! Add halo mask
-  afield = self%functionspaceInchalo%create_field(name='owned', kind=atlas_integer(kind(0)), levels=1)
-  allocate(hmask(self%isd:self%ied, self%jsd:self%jed))
-  hmask = 0
-  hmask(self%isc:self%iec, self%jsc:self%jec) = 1
-  call afield%data(int_ptr)
-  int_ptr(1,:) = pack(hmask, self%valid_halo_mask)
-  call afieldset%add(afield)
-  call afield%final()
+  ! halo exchanges for some of the fields 
+  call fGmask%halo_exchange()
+  call fInterpMask%halo_exchange()
+  call fArea%halo_exchange()
+  call fVertCoord%halo_exchange()
+  call fRossby%halo_exchange()
 
-  ! Add interpolation mask
-  afield = self%functionspaceInchalo%create_field(name='interp_mask', kind=atlas_real(kind_real), levels=1)
-  call afield%data(real_ptr)
-  real_ptr(1,:) = pack(self%mask2d, self%valid_halo_mask)
-  call afieldset%add(afield)
-  call afield%final()  
+  ! done, cleanup
+  call fArea%final()
+  call fInterpMask%final()
+  call fVertCoord%final()
+  call fGmask%final()
+  call fOwned%final()
+  call fRossby%final()
 
-end subroutine soca_geom_to_fieldset
+end subroutine soca_geom_init_fieldset
 
 
 ! ------------------------------------------------------------------------------
@@ -415,7 +381,7 @@ subroutine soca_geom_clone(self, other)
   self%rossby_radius = other%rossby_radius
   self%distance_from_coast = other%distance_from_coast
   self%h = other%h
-  self%valid_halo_mask = other%valid_halo_mask
+  self%atlas_ij2idx = other%atlas_ij2idx
   call self%fields_metadata%clone(other%fields_metadata)
 end subroutine soca_geom_clone
 
@@ -507,8 +473,6 @@ subroutine soca_geom_allocate(self)
   call soca_geom_get_domain_indices(self, "global", self%isg, self%ieg, self%jsg, self%jeg)
   call soca_geom_get_domain_indices(self, "compute", self%iscl, self%iecl, self%jscl, self%jecl, local=.true.)
   call soca_geom_get_domain_indices(self, "data", self%isdl, self%iedl, self%jsdl, self%jedl, local=.true.)
-  self%ngrid = (self%iec-self%isc+1) * (self%jec-self%jsc+1)
-  self%ngrid_halo = (self%ied-self%isd+1) * (self%jed-self%jsd+1)
 
   ! Allocate arrays on compute domain
   ! NOTE, sometimes MOM6 doesn't use all the halo points, we set 
@@ -538,8 +502,8 @@ subroutine soca_geom_allocate(self)
   allocate(self%rossby_radius(isd:ied,jsd:jed)); self%rossby_radius = 0.0_kind_real
   allocate(self%distance_from_coast(isd:ied,jsd:jed)); self%distance_from_coast = 0.0_kind_real
   allocate(self%h(isd:ied,jsd:jed,1:nzo));       self%h = 0.0_kind_real
-
-  allocate(self%valid_halo_mask(isd:ied,jsd:jed));self%valid_halo_mask = .true.
+  
+  allocate(self%atlas_ij2idx(isd:ied,jsd:jed));  self%atlas_ij2idx = -1
 
 end subroutine soca_geom_allocate
 
@@ -997,101 +961,105 @@ subroutine soca_geom_thickness2depth(self, h, z)
   end do
 end subroutine soca_geom_thickness2depth
 
-
 ! ------------------------------------------------------------------------------
-!> Copy a structured field into an ATLAS fieldset
-!!
-!! \related soca_geom_mod::soca_geom
-subroutine soca_geom_struct2atlas(self, dx_struct, dx_atlas)
-  class(soca_geom),     intent(in ) :: self
-  real(kind=kind_real), intent(in ) :: dx_struct(:,:)
-  type(fieldset_type),  intent(out) :: dx_atlas
-
-  real(kind_real), pointer :: real_ptr(:,:)
-  type(atlas_field) :: afield
-
-  dx_atlas = atlas_fieldset()
-  afield = self%functionspaceInchalo%create_field('var',kind=atlas_real(kind_real),levels=1)
-  call dx_atlas%add(afield)
-  call afield%data(real_ptr)
-  real_ptr(1,:) = pack(dx_struct, self%valid_halo_mask)
-  call afield%final()
-
-end subroutine soca_geom_struct2atlas
-
-
-! ------------------------------------------------------------------------------
-!> Copy a structured field from an ATLAS fieldset
-!!
-!! \related soca_geom_mod::soca_geom
-subroutine soca_geom_atlas2struct(self, dx_struct, dx_atlas)
-  class(soca_geom),     intent(in   ) :: self
-  real(kind=kind_real), intent(inout) :: dx_struct(:,:)
-  type(fieldset_type),  intent(inout) :: dx_atlas
-
-  real(kind_real), pointer :: real_ptr(:,:)
-  type(atlas_field) :: afield
-
-  afield = dx_atlas%field('var')
-  call afield%data(real_ptr)
-  dx_struct = unpack(real_ptr(1,:), self%valid_halo_mask, 0.0_kind_real)
-
-  call afield%final()
-
-end subroutine soca_geom_atlas2struct
-
-! ------------------------------------------------------------------------------
-!> Examine the halo gridpoints to determine which ones are non-duplicate and valid
-!!
-!! The OOPS interpolation does not like duplicate points, but MOM6 does have duplicate
-!! and/or unused halo points. We create a mask (self%valid_halo_mask) that is true if 1) 
-!! not a halo point, or 2) a halo point that is not a duplicate lat/lon already present 
-!! in the PE's halo or non-halo regions. Invalid halo points are also masked out 
-!! (i.e. the lat/lon are not updated from the initialed INVALID_HALO value when a halo 
-!! exchanged is performed
-subroutine soca_geom_find_invalid_halo(self)
-  class(soca_geom), intent(inout) :: self
-
-  integer ::  i1, i2, j1, j2
-
-  ! iterate over all halo gridpoints
-  outer: do j1 = self%jsd, self%jed
-    do i1 = self%isd, self%ied
-
-      ! if not halo region, skip
-      if (j1 >= self%jsc .and. j1 <= self%jec .and. &
-          i1 >= self%isc .and. i1 <= self%iec)  cycle
-
-      ! if gridpoint is still set to -999 after a halo exchange, it is an invalid point
-          if (self%lat(i1,j1) == INVALID_HALO) then
-            self%valid_halo_mask(i1,j1) = .false.
-            cycle
-          end if    
-          
-      ! make sure this halo point doesn't exist anywhere else, if it does mask it out
-      inner: do j2 = self%jsd, self%jed
-        do i2 = self%isd, self%ied
-          ! skip if point 1 is same as point 2
-          if (j1 == j2 .and. i1 == i2) cycle
-          ! skip if point 2 has already been masked out
-          if (.not. self%valid_halo_mask(i2,j2)) cycle
-          
-          ! if lat/lon are the same, mark point 1 as a duplicate, mask out
-          if (self%lat(i1,j1) == self%lat(i2,j2) .and. &
-              self%lon(i1,j1) == self%lon(i2,j2)) then
-            self%valid_halo_mask(i1,j1) = .false.
-            exit inner
-          end if
-        end do
-      end do inner
+! Get a 2d array of the valid nodes / cells that are to be used on this PE
+! for ATLAS mesh generation.
+! We assume that owned cells (quads) are generated north and east of the PE's owned nodes (vertices)
+! NOTE: this code assumes that a regional domain does NOT make use of halos on the boundary.
+!  We currently have our regional domain surrounded by "virtual land", so this is currently a correct
+!  assumption. When we get around to removing that "virtual land", and intend to populate the
+!  halo with a boundary conditions, the grid will need to be constructed slightly differently.
+subroutine soca_geom_mesh_valid_nodes_cells(self, nodes, cells)
+  class(soca_geom),   intent(inout) :: self
+  logical, allocatable, intent(out) :: nodes(:,:), cells(:,:)
+
+  integer :: i, j, start_i
+  logical :: tripolar, cyclic, regional
+
+  allocate(nodes(self%isc:self%iec+1, self%jsc:self%jec+1))
+  allocate(cells(self%isc:self%iec, self%jsc:self%jec))
+
+  nodes = .true.
+  cells = .true.
+
+  tripolar = iand(self%domain%Y_FLAGS, FOLD_NORTH_EDGE) /= 0
+  cyclic = iand(self%domain%X_FLAGS, CYCLIC_GLOBAL_DOMAIN) /= 0
+  regional = iand(self%domain%X_FLAGS, CYCLIC_GLOBAL_DOMAIN) == 0 .and. &
+             iand(self%domain%Y_FLAGS, CYCLIC_GLOBAL_DOMAIN) == 0
+
+  ! -------------------------------------------------------------------------------------
+  ! we need to mask out cells and/or nodes for 3 special cases.
+  ! Remove halo points that are otherwise going to be duplicated in the node list for the the PEs,
+  ! or, for regional (for now) remove unowned halo points at the domain boundary.
+  ! -------------------------------------------------------------------------------------
+
+  ! The grid is tripolar, and we are a PE at the northern edge
+  ! -------------------------------------------------------------------------------------
+  if(tripolar .and. self%jec == self%domain%NJGLOBAL) then
+      ! (start_i is the x index where we start removing nodes)
+    ! Remove all nodes in the eastern half.
+    start_i = max(self%domain%NIGLOBAL/2, self%isc)
+
+    ! additionally, this PE crosses the E/W midpoint, remove all extra nodes at the north for this PE
+    if (self%isc <= self%domain%NIGLOBAL/2 .and. self%iec > self%domain%NIGLOBAL/2) then      
+      start_i = self%isc
+    end if
     
+    ! remove some nodes
+    do i=start_i, self%iec+1
+        nodes(i, self%jec+1) = .false.
     end do
-  end do outer
 
-  ! count the number of points
-  self%ngrid_halo_valid = count(self%valid_halo_mask)  
+    ! remove some cells
+    do i=self%isc, self%iec
+      if (i >= self%domain%NIGLOBAL/2) then
+        cells(i, self%jec) = .false.
+      end if
+    end do
+  end if     
   
-end subroutine soca_geom_find_invalid_halo
+  ! the grid is cyclic in the E/W direction and we are a PE at the eastern edge
+  ! -------------------------------------------------------------------------------------  
+  if (cyclic .and. self%isc == 1 .and. self%iec == self%domain%NIGLOBAL) then  
+    ! cyclic boundaries, and this PE spans entire width.
+    ! Remove all nodes on the easternmost column.
+    do j=self%jsc, self%jec+1
+      nodes(self%iec+1, j) = .false.
+    end do
+  end if
+
+  ! We are a regional grid
+  ! NOTE: this logic should change at some point so that there ARE halos around the 
+  ! outer boundary
+  ! -------------------------------------------------------------------------------------
+  if (regional) then
+    ! are we are a PE at the top edge?
+    if(self%jec == self%domain%NJGLOBAL) then
+      ! remove nodes on top
+      do i=self%isc,self%iec+1
+        nodes(i, self%jec+1) = .false.
+      end do
+      ! remove cells on top
+      do i=self%isc, self%iec
+        cells(i, self%jec) = .false.
+      end do      
+    end if
+
+    ! are we a PE at the right edge?
+    if(self%iec == self%domain%NIGLOBAL) then
+      ! remove nodes on right
+      do j=self%jsc,self%jec+1
+        nodes(self%iec+1, j) = .false.
+      end do
+      ! remove cells on right
+      do j=self%jsc, self%jec
+        cells(self%iec, j) = .false.
+      end do          
+    end if
+
+  end if
+
+end subroutine
 
 ! ------------------------------------------------------------------------------
 
diff --git a/src/soca/Increment/Increment.cc b/src/soca/Increment/Increment.cc
index 895cc4f1d..ad501d492 100755
--- a/src/soca/Increment/Increment.cc
+++ b/src/soca/Increment/Increment.cc
@@ -1,5 +1,5 @@
 /*
- * (C) Copyright 2017-2022 UCAR
+ * (C) Copyright 2017-2023 UCAR
  *
  * This software is licensed under the terms of the Apache Licence Version 2.0
  * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
@@ -358,13 +358,7 @@ namespace soca {
 
   void Increment::toFieldSet(atlas::FieldSet &fs) const {
     soca_increment_to_fieldset_f90(toFortran(), vars_, fs.get());
-  }
-
-// -----------------------------------------------------------------------------
-
-  void Increment::toFieldSetAD(const atlas::FieldSet &fs) {
-    if (fs.empty()) return;
-    soca_increment_to_fieldset_ad_f90(toFortran(), vars_, fs.get());
+    geom_.functionSpace().haloExchange(fs);
   }
 
 // -----------------------------------------------------------------------------
diff --git a/src/soca/Increment/Increment.h b/src/soca/Increment/Increment.h
index 4c4df95df..0d25add8a 100644
--- a/src/soca/Increment/Increment.h
+++ b/src/soca/Increment/Increment.h
@@ -1,6 +1,6 @@
 /*
  * (C) Copyright 2009-2016 ECMWF.
- * (C) Copyright 2017-2022 UCAR.
+ * (C) Copyright 2017-2023 UCAR.
  *
  * This software is licensed under the terms of the Apache Licence Version 2.0
  * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
@@ -90,7 +90,6 @@ namespace soca {
 
       /// ATLAS
       void toFieldSet(atlas::FieldSet &) const;
-      void toFieldSetAD(const atlas::FieldSet &);
       void fromFieldSet(const atlas::FieldSet &);
 
       /// I/O and diagnostics
diff --git a/src/soca/Increment/IncrementFortran.h b/src/soca/Increment/IncrementFortran.h
index 6b09e5cce..2bc722ddd 100644
--- a/src/soca/Increment/IncrementFortran.h
+++ b/src/soca/Increment/IncrementFortran.h
@@ -1,5 +1,5 @@
 /*
- * (C) Copyright 2020-2022 UCAR
+ * (C) Copyright 2020-2023 UCAR
  *
  * This software is licensed under the terms of the Apache Licence Version 2.0
  * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
@@ -60,9 +60,6 @@ namespace soca {
     void soca_increment_to_fieldset_f90(const F90flds &,
                                         const oops::Variables &,
                                         atlas::field::FieldSetImpl *);
-    void soca_increment_to_fieldset_ad_f90(const F90flds &,
-                                           const oops::Variables &,
-                                           const atlas::field::FieldSetImpl *);
     void soca_increment_from_fieldset_f90(const F90flds &,
                                    const oops::Variables &,
                                    const atlas::field::FieldSetImpl *);
diff --git a/src/soca/Increment/soca_increment.interface.F90 b/src/soca/Increment/soca_increment.interface.F90
index 1112d7500..310d8c531 100644
--- a/src/soca/Increment/soca_increment.interface.F90
+++ b/src/soca/Increment/soca_increment.interface.F90
@@ -1,4 +1,4 @@
-! (C) Copyright 2020-2022 UCAR
+! (C) Copyright 2020-2023 UCAR
 !
 ! This software is licensed under the terms of the Apache Licence Version 2.0
 ! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
@@ -633,24 +633,5 @@ subroutine soca_increment_vert_scales_c(c_key_self, c_vert) bind(c,name='soca_in
 
 end subroutine soca_increment_vert_scales_c
 
-! ------------------------------------------------------------------------------
-!> C++ interface for Increment version of soca_field_mod::soca_field::get_fieldset_ad()
-subroutine soca_increment_to_fieldset_ad_c(c_key_self, c_vars, c_fieldset) &
-    bind (c, name='soca_increment_to_fieldset_ad_f90')
-  integer(c_int),       intent(in) :: c_key_self
-  type(c_ptr), value,   intent(in) :: c_vars
-  type(c_ptr), value,   intent(in) :: c_fieldset
-
-  type(soca_increment), pointer :: self
-  type(oops_variables) :: vars
-  type(atlas_fieldset) :: afieldset
-
-  call soca_increment_registry%get(c_key_self, self)
-  vars = oops_variables(c_vars)
-  afieldset = atlas_fieldset(c_fieldset)
-
-  call self%to_fieldset_ad(vars, afieldset)
-end subroutine
-
 ! ------------------------------------------------------------------------------
 end module
diff --git a/src/soca/LinearVariableChange/HorizFilt/soca_horizfilt_mod.F90 b/src/soca/LinearVariableChange/HorizFilt/soca_horizfilt_mod.F90
index 5a9b1dac6..b8c176799 100644
--- a/src/soca/LinearVariableChange/HorizFilt/soca_horizfilt_mod.F90
+++ b/src/soca/LinearVariableChange/HorizFilt/soca_horizfilt_mod.F90
@@ -1,4 +1,4 @@
-! (C) Copyright 2017-2021 UCAR.
+! (C) Copyright 2017-2023 UCAR.
 !
 ! This software is licensed under the terms of the Apache Licence Version 2.0
 ! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
@@ -109,11 +109,6 @@ subroutine soca_horizfilt_setup(self, f_conf, geom, traj, vars)
       do i = self%isc, self%iec
         do ii = -1,1
             do jj = -1,1
-              ! if the adjacent cell is an invalid halo, skip it
-              if (.not. geom%valid_halo_mask(i+ii, j+jj)) then
-                dist(ii, jj) = 0.0
-                cycle
-              end if
 
               ! Great circle distance
               if(self%scale_dist > 0) then
diff --git a/src/soca/State/State.cc b/src/soca/State/State.cc
index 49e4ce931..1ef2b1af8 100644
--- a/src/soca/State/State.cc
+++ b/src/soca/State/State.cc
@@ -1,5 +1,5 @@
 /*
- * (C) Copyright 2017-2021 UCAR
+ * (C) Copyright 2017-2023 UCAR
  *
  * This software is licensed under the terms of the Apache Licence Version 2.0
  * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
@@ -255,6 +255,7 @@ namespace soca {
 
   void State::toFieldSet(atlas::FieldSet &fset) const {
     soca_state_to_fieldset_f90(toFortran(), vars_, fset.get());
+    fset.haloExchange();
   }
 
   // -----------------------------------------------------------------------------
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 0beebbd20..0c10cf7d8 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -613,8 +613,10 @@ soca_add_test( NAME dirac_soca_nomask
                TEST_DEPENDS test_soca_gridgen
                             test_soca_static_socaerror_init)
 
+# NOTE the following test throws a FPE only on clang
 soca_add_test( NAME dirac_diffusion
                EXE soca_error_covariance_toolbox.x
+               NOTRAPFPE
                TEST_DEPENDS test_soca_parameters_diffusion)
 
 # h(x) executables
diff --git a/test/testinput/gridgen.yml b/test/testinput/gridgen.yml
index cc5769f7f..80fc88575 100644
--- a/test/testinput/gridgen.yml
+++ b/test/testinput/gridgen.yml
@@ -6,6 +6,10 @@ geometry:
   fields metadata: data_static/fields_metadata.yml
   rossby file: data_static/rossrad.dat
 
+  # save the atlas mesh for viewing in gmsh
+  gmsh save: true
+  gmsh filename: data_generated/gridgen/gmsh.msh
+
 test:
   reference filename: testref/gridgen.test
   test output filename: testoutput/gridgen.test
diff --git a/test/testref/3dhyb.test b/test/testref/3dhyb.test
index ed031e925..fe833233d 100644
--- a/test/testref/3dhyb.test
+++ b/test/testref/3dhyb.test
@@ -1,17 +1,17 @@
 CostJb   : Nonlinear Jb = 0.0000000000000000e+00
-CostJo   : Nonlinear Jo(ADT) = 2.0458137413924788e+02, nobs = 99, Jo/n = 2.0664785266590697e+00, err = 1.0000000149011613e-01
-CostFunction: Nonlinear J = 2.0458137413924788e+02
-RPCGMinimizer: reduction in residual norm = 2.1404205571904238e-02
+CostJo   : Nonlinear Jo(ADT) = 2.0980804284155522e+02, nobs = 99, Jo/n = 2.1192731600157093e+00, err = 1.0000000149011613e-01
+CostFunction: Nonlinear J = 2.0980804284155522e+02
+RPCGMinimizer: reduction in residual norm = 1.3986916105041314e-02
 CostFunction::addIncrement: Analysis: 
   Valid time: 2018-04-15T00:00:00Z
-  cicen   min=-0.0040988630419505   max=1.0000630744489174   mean=0.1175050968515152
-  hicen   min=-0.1725638320233036   max=4.0326711168772373   mean=0.4671340927388408
+  cicen   min=-0.0055833351906180   max=1.0000625277294939   mean=0.1175137930143875
+  hicen   min=-0.1340300448550523   max=4.0326711192074995   mean=0.4668708457595877
   hsnon   min=0.0000000000000000   max=1.2712833951042413   mean=0.0886865877527975
-   socn   min=10.7210460395083924   max=40.4416591897031168   mean=34.5753803016227792
-   tocn   min=-1.8973132653785527   max=31.7004647372264117   mean=6.0143150203319413
-   uocn   min=-0.8581810504448349   max=0.7001468341842769   mean=-0.0002551821852330
-   vocn   min=-0.7656654371938378   max=1.4377766424239879   mean=0.0021996485608986
-    ssh   min=-2.1639526837757632   max=0.8962563828063304   mean=-0.2944013250972501
+   socn   min=10.7210460395083924   max=40.4416591897031168   mean=34.5771060899726166
+   tocn   min=-1.8977624967437274   max=31.7004647317691663   mean=6.0154772032339441
+   uocn   min=-0.8581811478246130   max=0.7001477830089919   mean=-0.0002559066432762
+   vocn   min=-0.7656861309282604   max=1.4377766510613166   mean=0.0021994321716903
+    ssh   min=-2.1760245287885742   max=0.8999953604193408   mean=-0.2948395432550496
    hocn   min=0.0009999999999977   max=1345.6400000000003274   mean=128.6280642065023017
      sw   min=-225.0977630615234375   max=-0.0000000000000000   mean=-71.7393205198725354
     lhf   min=-38.1373977661132812   max=256.5976562500000000   mean=42.0100402226406899
@@ -23,6 +23,6 @@ layer_depth   min=2.2854716757984130   max=5658.3057467114012979   mean=1200.522
 
 
 
-CostJb   : Nonlinear Jb = 6.0654922897128243
-CostJo   : Nonlinear Jo(ADT) = 115.5460048800289599, nobs = 99, Jo/n = 1.1671313624245350, err = 0.1000000014901161
-CostFunction: Nonlinear J = 121.6114971697417815
+CostJb   : Nonlinear Jb = 6.2153370309095379
+CostJo   : Nonlinear Jo(ADT) = 117.2133340087334972, nobs = 99, Jo/n = 1.1839730707952878, err = 0.1000000014901161
+CostFunction: Nonlinear J = 123.4286710396430351
diff --git a/test/testref/3dhybfgat.test b/test/testref/3dhybfgat.test
index 43f8b968f..e1a492318 100644
--- a/test/testref/3dhybfgat.test
+++ b/test/testref/3dhybfgat.test
@@ -1,16 +1,16 @@
 CostJb   : Nonlinear Jb = 0.0000000000000000e+00
-CostJo   : Nonlinear Jo(ADT) = 8.7370280580244312e+01, nobs = 31, Jo/n = 2.8183961477498167e+00, err = 1.0000000149011612e-01
-CostFunction: Nonlinear J = 8.7370280580244312e+01
-RPCGMinimizer: reduction in residual norm = 3.0018932255325439e-03
+CostJo   : Nonlinear Jo(ADT) = 8.7827980335709171e+01, nobs = 31, Jo/n = 2.8331606559906186e+00, err = 1.0000000149011612e-01
+CostFunction: Nonlinear J = 8.7827980335709171e+01
+RPCGMinimizer: reduction in residual norm = 2.2743017181226293e-03
 CostFunction::addIncrement: Analysis: 
   Valid time: 2018-04-15T00:00:00Z
-  cicen   min=-0.0040535217089970   max=1.0000000421973050   mean=0.1174810780852600
-  hicen   min=-0.0902428957660459   max=4.0326673084246947   mean=0.4697245454894248
-   socn   min=10.7210460395083924   max=40.4416591897031168   mean=34.5590151339643867
-   uocn   min=-0.8581694001038450   max=0.7001584930307920   mean=-0.0002536521284757
-   vocn   min=-0.7661101225938430   max=1.4377766409421606   mean=0.0021979270760311
-   tocn   min=-1.8883899367449197   max=31.7004645720658580   mean=6.0142205656355632
-    ssh   min=-1.9233869940455874   max=0.8888169925302998   mean=-0.2881240363183634
+  cicen   min=-0.0055094416233540   max=1.0000000420883628   mean=0.1174900878620023
+  hicen   min=-0.0898682624110847   max=4.0326673084246947   mean=0.4695679968317440
+   socn   min=10.7210460395083924   max=40.4416591897031168   mean=34.5587872384660102
+   uocn   min=-0.8581694003631951   max=0.7001563851954871   mean=-0.0002538829121215
+   vocn   min=-0.7661101020702997   max=1.4377766409421606   mean=0.0021980693625999
+   tocn   min=-1.8883899367432375   max=31.7004645720658580   mean=6.0145488739534851
+    ssh   min=-1.9232987964448847   max=0.8928791626958916   mean=-0.2875104304210939
    hocn   min=0.0009999999999977   max=1345.6400000000003274   mean=128.6280642065023017
      sw   min=0.0000000000000000   max=0.0000000000000000   mean=0.0000000000000000
     lhf   min=0.0000000000000000   max=0.0000000000000000   mean=0.0000000000000000
@@ -22,6 +22,6 @@ layer_depth   min=2.2854716757984130   max=5658.3057467114012979   mean=1200.522
 
 
 
-CostJb   : Nonlinear Jb = 2.2768042398308497
-CostJo   : Nonlinear Jo(ADT) = 77.2923955680257535, nobs = 31, Jo/n = 2.4933030828395406, err = 0.1000000014901161
-CostFunction: Nonlinear J = 79.5691998078566058
+CostJb   : Nonlinear Jb = 2.2400054920788746
+CostJo   : Nonlinear Jo(ADT) = 76.5887027511879666, nobs = 31, Jo/n = 2.4706033145544506, err = 0.1000000014901161
+CostFunction: Nonlinear J = 78.8287082432668456
diff --git a/test/testref/3dvar_godas.test b/test/testref/3dvar_godas.test
index 8cc99ee13..34e656bec 100644
--- a/test/testref/3dvar_godas.test
+++ b/test/testref/3dvar_godas.test
@@ -1,27 +1,27 @@
 CostJb   : Nonlinear Jb = 0.0000000000000000e+00
-CostJo   : Nonlinear Jo(CoolSkin) = 1.6731214283333800e+04, nobs = 200, Jo/n = 8.3656071416669008e+01, err = 3.6419991997020579e-01
-CostJo   : Nonlinear Jo(SeaSurfaceTemp) = 1.6258116856443858e+03, nobs = 174, Jo/n = 9.3437453197953211e+00, err = 3.6280357515271239e-01
-CostJo   : Nonlinear Jo(SeaSurfaceSalinity) = 5.4431675889674533e+00, nobs = 29, Jo/n = 1.8769543410232598e-01, err = 1.0000000000000000e+00
-CostJo   : Nonlinear Jo(ADT) = 1.0375372686131383e-02, nobs = 95, Jo/n = 1.0921444932769877e-04, err = 2.8169689856196570e+01
-CostJo   : Nonlinear Jo(InsituTemperature) = 3.3944990007872309e+02, nobs = 205, Jo/n = 1.6558531711157225e+00, err = 8.9723213032395277e-01
-CostJo   : Nonlinear Jo(InsituSalinity) = 3.4804589303321143e+02, nobs = 218, Jo/n = 1.5965407937303278e+00, err = 6.0819699545777528e-01
-CostJo   : Nonlinear Jo(SeaIceFraction) = 5.8152359560103514e+02, nobs = 88, Jo/n = 6.6082226772844903e+00, err = 1.0000000149011612e-01
-CostFunction: Nonlinear J = 1.9631498900652809e+04
-RPCGMinimizer: reduction in residual norm = 1.4894433280145411e-01
+CostJo   : Nonlinear Jo(CoolSkin) = 1.6329134149677742e+04, nobs = 199, Jo/n = 8.2055950500893175e+01, err = 3.6419991997020579e-01
+CostJo   : Nonlinear Jo(SeaSurfaceTemp) = 1.6066024929416167e+03, nobs = 173, Jo/n = 9.2867196123792866e+00, err = 3.6322734119159583e-01
+CostJo   : Nonlinear Jo(SeaSurfaceSalinity) = 5.5136316271599295e+00, nobs = 29, Jo/n = 1.9012522852275618e-01, err = 1.0000000000000000e+00
+CostJo   : Nonlinear Jo(ADT) = 1.0334301354601411e-02, nobs = 95, Jo/n = 1.0878211952212012e-04, err = 2.8230815589406479e+01
+CostJo   : Nonlinear Jo(InsituTemperature) = 3.4020248029963227e+02, nobs = 205, Jo/n = 1.6595242941445476e+00, err = 8.9723213032395277e-01
+CostJo   : Nonlinear Jo(InsituSalinity) = 3.4727192955958594e+02, nobs = 218, Jo/n = 1.5929905025669080e+00, err = 6.0819699545777528e-01
+CostJo   : Nonlinear Jo(SeaIceFraction) = 5.6944306282109505e+02, nobs = 89, Jo/n = 6.3982366609111807e+00, err = 1.0000000149011612e-01
+CostFunction: Nonlinear J = 1.9198178081228183e+04
+RPCGMinimizer: reduction in residual norm = 1.3394869278655661e-01
 CostFunction::addIncrement: Analysis: 
   Valid time: 2018-04-15T00:00:00Z
-  cicen   min=-0.0000705792064587   max=1.0003375268744732   mean=0.1175473153867414
+  cicen   min=-0.0000563665143372   max=1.0003307403191071   mean=0.1175499389485267
   hicen   min=0.0000000000000000   max=4.0326673084246947   mean=0.4712515705773916
   hsnon   min=0.0000000000000000   max=1.2712833951042413   mean=0.0886865877527975
-   socn   min=10.7210460395083924   max=40.4416591897031168   mean=34.5443924231097128
-   tocn   min=-1.8883899372702533   max=31.7004645720658580   mean=6.0183795308274490
-    ssh   min=-1.9244192301102359   max=0.9272827055997203   mean=-0.2767862045255065
+   socn   min=10.7210460395083924   max=40.4416591897031168   mean=34.5443925948116970
+   tocn   min=-1.8883899372702533   max=31.7004645720658580   mean=6.0183977449112982
+    ssh   min=-1.9244089651560510   max=0.9272827138726495   mean=-0.2767862613297231
    hocn   min=0.0009999999999977   max=1345.6400000000003274   mean=128.6280642065023017
-     sw   min=-225.0942669976566890   max=0.0000000000000000   mean=-71.7383324463225733
-    lhf   min=-38.1375347797083535   max=256.5954640460074643   mean=42.0097100740689484
-    shf   min=-58.9492919026422513   max=201.3742668398652995   mean=6.0753281714416643
-     lw   min=0.0000000000000000   max=88.0360937135157400   mean=31.3955347883638041
-     us   min=0.0044037898505438   max=0.0196682633764925   mean=0.0086888474336386
+     sw   min=-225.0943644800591414   max=0.0000000000000000   mean=-71.7382903028821204
+    lhf   min=-38.1375465556171065   max=256.5949994230312541   mean=42.0096946792278629
+    shf   min=-58.9492928606160191   max=201.3742494565976813   mean=6.0753275791560339
+     lw   min=0.0000000000000000   max=88.0360940812700505   mean=31.3955349934081234
+     us   min=0.0044042662455111   max=0.0196689588586444   mean=0.0086900813793515
     mld   min=2.2854716757984130   max=4593.1533423819937525   mean=192.4109073940401515
 layer_depth   min=2.2854716757984130   max=5658.3057467114012979   mean=1200.5229536158392420
 
@@ -33,12 +33,12 @@ layer_depth   min=2.2854716757984130   max=5658.3057467114012979   mean=1200.522
 
 
 
-CostJb   : Nonlinear Jb = 0.4007664673415510
-CostJo   : Nonlinear Jo(CoolSkin) = 16185.1504262224261765, nobs = 200, Jo/n = 80.9257521311121337, err = 0.3641999199702058
-CostJo   : Nonlinear Jo(SeaSurfaceTemp) = 1596.8311924519289278, nobs = 174, Jo/n = 9.1771907612179824, err = 0.3628035751527124
-CostJo   : Nonlinear Jo(SeaSurfaceSalinity) = 5.4426955883488013, nobs = 29, Jo/n = 0.1876791582189242, err = 1.0000000000000000
-CostJo   : Nonlinear Jo(ADT) = 0.0103757275797161, nobs = 95, Jo/n = 0.0001092181850496, err = 28.1696898561965696
-CostJo   : Nonlinear Jo(InsituTemperature) = 337.3241194423258094, nobs = 205, Jo/n = 1.6454835094747600, err = 0.8972321303239528
-CostJo   : Nonlinear Jo(InsituSalinity) = 346.7332613876146752, nobs = 218, Jo/n = 1.5905195476496086, err = 0.6081969954577753
-CostJo   : Nonlinear Jo(SeaIceFraction) = 580.7063304346302175, nobs = 88, Jo/n = 6.5989355731207979, err = 0.1000000014901161
-CostFunction: Nonlinear J = 19052.5991677221936698
+CostJb   : Nonlinear Jb = 0.4158118826125912
+CostJo   : Nonlinear Jo(CoolSkin) = 15785.5145342027608422, nobs = 199, Jo/n = 79.3241936392098523, err = 0.3641999199702058
+CostJo   : Nonlinear Jo(SeaSurfaceTemp) = 1574.3766446106681087, nobs = 173, Jo/n = 9.1004430324316079, err = 0.3632273411915958
+CostJo   : Nonlinear Jo(SeaSurfaceSalinity) = 5.5131254646506047, nobs = 29, Jo/n = 0.1901077746431243, err = 1.0000000000000000
+CostJo   : Nonlinear Jo(ADT) = 0.0103346677169496, nobs = 95, Jo/n = 0.0001087859759679, err = 28.2308155894064789
+CostJo   : Nonlinear Jo(InsituTemperature) = 337.5858413149268245, nobs = 205, Jo/n = 1.6467602015362284, err = 0.8972321303239528
+CostJo   : Nonlinear Jo(InsituSalinity) = 345.8714449969637599, nobs = 218, Jo/n = 1.5865662614539622, err = 0.6081969954577753
+CostJo   : Nonlinear Jo(SeaIceFraction) = 568.6042208009064325, nobs = 89, Jo/n = 6.3888114696731062, err = 0.1000000014901161
+CostFunction: Nonlinear J = 18617.8919579412031453
diff --git a/test/testref/3dvar_soca.test b/test/testref/3dvar_soca.test
index 54e3816a8..d33000762 100644
--- a/test/testref/3dvar_soca.test
+++ b/test/testref/3dvar_soca.test
@@ -1,31 +1,31 @@
 CostJb   : Nonlinear Jb = 0.0000000000000000e+00
-CostJo   : Nonlinear Jo(CoolSkin) = 1.4169933330087919e+04, nobs = 174, Jo/n = 8.1436398448781148e+01, err = 3.6280357515271239e-01
-CostJo   : Nonlinear Jo(SeaSurfaceTemp) = 7.7388800526767272e+02, nobs = 147, Jo/n = 5.2645442535215832e+00, err = 3.6246510156485134e-01
-CostJo   : Nonlinear Jo(SeaSurfaceSalinity) = 5.4431675889674533e+00, nobs = 29, Jo/n = 1.8769543410232598e-01, err = 1.0000000000000000e+00
-CostJo   : Nonlinear Jo(ADT) = 1.3650926098993568e+02, nobs = 91, Jo/n = 1.5001017691201723e+00, err = 1.0000000149011612e-01
-CostJo   : Nonlinear Jo(InsituTS) = 6.8749579311193474e+02, nobs = 423, Jo/n = 1.6252855629123752e+00, err = 7.6208809857781845e-01
-CostJo   : Nonlinear Jo(SeaIceFraction) = 5.8152359560103514e+02, nobs = 88, Jo/n = 6.6082226772844903e+00, err = 1.0000000149011612e-01
-CostJo   : Nonlinear Jo(SeaSurfaceChlorophyll) = 1.4791913800716829e+03, nobs = 129, Jo/n = 1.1466599845516923e+01, err = 6.1865027747001339e-02
+CostJo   : Nonlinear Jo(CoolSkin) = 1.4174299653183771e+04, nobs = 173, Jo/n = 8.1932367937478446e+01, err = 3.6322734119159583e-01
+CostJo   : Nonlinear Jo(SeaSurfaceTemp) = 7.7377196587412891e+02, nobs = 146, Jo/n = 5.2998079854392390e+00, err = 3.6296533314631141e-01
+CostJo   : Nonlinear Jo(SeaSurfaceSalinity) = 5.5136316271599295e+00, nobs = 29, Jo/n = 1.9012522852275618e-01, err = 1.0000000000000000e+00
+CostJo   : Nonlinear Jo(ADT) = 1.3879704103180882e+02, nobs = 91, Jo/n = 1.5252422091407563e+00, err = 1.0000000149011612e-01
+CostJo   : Nonlinear Jo(InsituTS) = 6.8747440985921844e+02, nobs = 423, Jo/n = 1.6252350114875140e+00, err = 7.6208809857781845e-01
+CostJo   : Nonlinear Jo(SeaIceFraction) = 5.6944306282109505e+02, nobs = 89, Jo/n = 6.3982366609111807e+00, err = 1.0000000149011612e-01
+CostJo   : Nonlinear Jo(SeaSurfaceChlorophyll) = 1.3639731179097098e+03, nobs = 129, Jo/n = 1.0573435022555890e+01, err = 6.1865027747001339e-02
 CostJo   : Nonlinear Jo(SeaSurfaceBiomassP) = 4.1314922104098372e+02, nobs = 112, Jo/n = 3.6888323307230690e+00, err = 1.1692603016453901e-08
 CostJo   : Nonlinear Jo(SeaSurfaceHS) = 1.0896018677668069e+03, nobs = 6, Jo/n = 1.8160031129446782e+02, err = 1.0000000149011612e-01
-CostFunction: Nonlinear J = 1.9336735621526943e+04
-RPCGMinimizer: reduction in residual norm = 1.4943903814260545e-01
+CostFunction: Nonlinear J = 1.9216023971114686e+04
+RPCGMinimizer: reduction in residual norm = 1.3357822017190246e-01
 CostFunction::addIncrement: Analysis: 
   Valid time: 2018-04-15T00:00:00Z
-  cicen   min=-0.0001983484661628   max=1.0004563903533681   mean=0.1175584987640446
+  cicen   min=-0.0002115179073194   max=1.0004820341037177   mean=0.1175659422442119
   hicen   min=0.0000000000000000   max=4.0326673084246947   mean=0.4712515705773916
-   socn   min=10.7210460395083924   max=40.4416591897031168   mean=34.5444154414926672
-   tocn   min=-1.8883899372702533   max=31.7004645720658580   mean=6.0193560753976145
-    ssh   min=-1.9246079806628480   max=0.9272948016729473   mean=-0.2767213847127198
+   socn   min=10.7210460395083924   max=40.4416591897031168   mean=34.5444197182825405
+   tocn   min=-1.8883899372702533   max=31.7004645720658580   mean=6.0194012222162385
+    ssh   min=-1.9246315345115299   max=0.9272891664663706   mean=-0.2767204306215251
    hocn   min=0.0009999999999977   max=1345.6400000000003274   mean=128.6280642065023017
-     sw   min=-225.0945898286534543   max=0.0000000000000000   mean=-71.7382190891333664
-    lhf   min=-38.1375526767112092   max=256.5953230482455183   mean=42.0096679932438946
-    shf   min=-58.9492904310832628   max=201.3742561297929115   mean=6.0753266889941919
-     lw   min=0.0000000000000000   max=88.0360960036400826   mean=31.3955354121763648
-     us   min=0.0044062523598695   max=0.0196718967384776   mean=0.0086921922934188
-    chl   min=-0.0000132860100953   max=4.6719899196170607   mean=0.1184787259407491
-   biop   min=0.0000000000000000   max=0.0000001868533390   mean=0.0000000065395151
-    swh   min=0.0000996535891318   max=6.7491955757141104   mean=1.0621519653888534
+     sw   min=-225.0946147629603331   max=0.0000000000000000   mean=-71.7381266251842646
+    lhf   min=-38.1375723103176085   max=256.5947019639968403   mean=42.0096297492159820
+    shf   min=-58.9492914708947495   max=201.3742276461881602   mean=6.0753253636382611
+     lw   min=0.0000000000000000   max=88.0360972107224740   mean=31.3955360261616576
+     us   min=0.0044072290018029   max=0.0196739899437386   mean=0.0086954229389514
+    chl   min=-0.0000081536937266   max=4.6719899196170607   mean=0.1184780720546971
+   biop   min=0.0000000000000000   max=0.0000001868533390   mean=0.0000000065395206
+    swh   min=0.0000996535891318   max=6.7491955757141104   mean=1.0621545944148620
     mld   min=2.2854716757984130   max=4593.1533423819937525   mean=192.4109073940401515
 layer_depth   min=2.2854716757984130   max=5658.3057467114012979   mean=1200.5229536158392420
 
@@ -39,14 +39,14 @@ layer_depth   min=2.2854716757984130   max=5658.3057467114012979   mean=1200.522
 
 
 
-CostJb   : Nonlinear Jb = 3.9507913202249219
-CostJo   : Nonlinear Jo(CoolSkin) = 13655.1163728023566364, nobs = 174, Jo/n = 78.4776803034618240, err = 0.3628035751527124
-CostJo   : Nonlinear Jo(SeaSurfaceTemp) = 769.1250325174885347, nobs = 147, Jo/n = 5.2321430783502620, err = 0.3624651015648513
-CostJo   : Nonlinear Jo(SeaSurfaceSalinity) = 5.4427739531973049, nobs = 29, Jo/n = 0.1876818604550795, err = 1.0000000000000000
-CostJo   : Nonlinear Jo(ADT) = 136.4788599610268420, nobs = 91, Jo/n = 1.4997676918794158, err = 0.1000000014901161
-CostJo   : Nonlinear Jo(InsituTS) = 686.2447930167938921, nobs = 423, Jo/n = 1.6223281158789453, err = 0.7620880985778185
-CostJo   : Nonlinear Jo(SeaIceFraction) = 580.4174851187933655, nobs = 88, Jo/n = 6.5956532399862882, err = 0.1000000014901161
-CostJo   : Nonlinear Jo(SeaSurfaceChlorophyll) = 1460.0295525378164712, nobs = 129, Jo/n = 11.3180585468047781, err = 0.0618650277470013
-CostJo   : Nonlinear Jo(SeaSurfaceBiomassP) = 409.2402710143757645, nobs = 112, Jo/n = 3.6539309911997835, err = 0.0000000116926030
-CostJo   : Nonlinear Jo(SeaSurfaceHS) = 1086.9456348424896532, nobs = 6, Jo/n = 181.1576058070816089, err = 0.1000000014901161
-CostFunction: Nonlinear J = 18792.9915670845621207
+CostJb   : Nonlinear Jb = 4.2941089565347577
+CostJo   : Nonlinear Jo(CoolSkin) = 13638.9679733031734941, nobs = 173, Jo/n = 78.8379651636021634, err = 0.3632273411915958
+CostJo   : Nonlinear Jo(SeaSurfaceTemp) = 768.2386717644395731, nobs = 146, Jo/n = 5.2619087107153399, err = 0.3629653331463114
+CostJo   : Nonlinear Jo(SeaSurfaceSalinity) = 5.5131762689470261, nobs = 29, Jo/n = 0.1901095265154147, err = 1.0000000000000000
+CostJo   : Nonlinear Jo(ADT) = 138.7634066289171528, nobs = 91, Jo/n = 1.5248726003177708, err = 0.1000000014901161
+CostJo   : Nonlinear Jo(InsituTS) = 686.0324173035821786, nobs = 423, Jo/n = 1.6218260456349460, err = 0.7620880985778185
+CostJo   : Nonlinear Jo(SeaIceFraction) = 568.2192882044269027, nobs = 89, Jo/n = 6.3844863843194037, err = 0.1000000014901161
+CostJo   : Nonlinear Jo(SeaSurfaceChlorophyll) = 1344.5348295591829810, nobs = 129, Jo/n = 10.4227506167378525, err = 0.0618650277470013
+CostJo   : Nonlinear Jo(SeaSurfaceBiomassP) = 408.6554750400586045, nobs = 112, Jo/n = 3.6487095985719518, err = 0.0000000116926030
+CostJo   : Nonlinear Jo(SeaSurfaceHS) = 1086.5431097938458151, nobs = 6, Jo/n = 181.0905182989743025, err = 0.1000000014901161
+CostFunction: Nonlinear J = 18649.7624568231112789
diff --git a/test/testref/3dvarbump.test b/test/testref/3dvarbump.test
index beb9ce12f..e00afe83b 100644
--- a/test/testref/3dvarbump.test
+++ b/test/testref/3dvarbump.test
@@ -1,19 +1,19 @@
 CostJb   : Nonlinear Jb = 0.0000000000000000e+00
-CostJo   : Nonlinear Jo(SeaSurfaceTemp) = 7.2761520714484413e+02, nobs = 123, Jo/n = 5.9155707897954803e+00, err = 3.7062899283957257e-01
-CostJo   : Nonlinear Jo(SeaSurfaceSalinity) = 5.7834714574083208e+01, nobs = 50, Jo/n = 1.1566942914816642e+00, err = 1.0000000000000000e+00
-CostJo   : Nonlinear Jo(ADT) = 1.7801726582342692e+02, nobs = 91, Jo/n = 1.9562336903673287e+00, err = 1.0000000149011612e-01
-CostJo   : Nonlinear Jo(InsituTemperature) = 3.4603234526181200e+02, nobs = 203, Jo/n = 1.7045928338020295e+00, err = 9.0034813852418905e-01
-CostJo   : Nonlinear Jo(InsituSalinity) = 3.1204797269182507e+02, nobs = 218, Jo/n = 1.4314127187698398e+00, err = 6.0819699545777528e-01
-CostFunction: Nonlinear J = 1.6215475054959916e+03
-RPCGMinimizer: reduction in residual norm = 1.8327341126136407e-01
+CostJo   : Nonlinear Jo(SeaSurfaceTemp) = 7.1091386479968503e+02, nobs = 123, Jo/n = 5.7797875186966259e+00, err = 3.7098966299006642e-01
+CostJo   : Nonlinear Jo(SeaSurfaceSalinity) = 5.8283270347772934e+01, nobs = 50, Jo/n = 1.1656654069554586e+00, err = 1.0000000000000000e+00
+CostJo   : Nonlinear Jo(ADT) = 1.9163642890807802e+02, nobs = 93, Jo/n = 2.0606067624524518e+00, err = 1.0000000149011612e-01
+CostJo   : Nonlinear Jo(InsituTemperature) = 3.9116568630186021e+02, nobs = 205, Jo/n = 1.9081252990334645e+00, err = 9.0169108266903031e-01
+CostJo   : Nonlinear Jo(InsituSalinity) = 3.1078118440821123e+02, nobs = 218, Jo/n = 1.4256017633404185e+00, err = 6.0819699545777528e-01
+CostFunction: Nonlinear J = 1.6627804347656074e+03
+RPCGMinimizer: reduction in residual norm = 1.6180170436945335e-01
 CostFunction::addIncrement: Analysis: 
   Valid time: 2018-04-15T00:00:00Z
    hocn   min=0.0009999999999977   max=1345.6400000000003274   mean=132.2975966204909355
-   socn   min=21.6778874359345082   max=39.8774684135234310   mean=34.6314659731017684
-   tocn   min=-1.8883899372702533   max=32.9461357419237331   mean=6.3373244093251282
+   socn   min=21.6778874359345082   max=39.8774684135234310   mean=34.6333415889160463
+   tocn   min=-1.9147364365672823   max=32.5667623078130362   mean=6.3366284689423633
    uocn   min=-0.6975788196638589   max=0.4840625031242498   mean=0.0013055310445162
    vocn   min=-0.7661101215480253   max=0.9402890903466523   mean=0.0028998456865158
-    ssh   min=-2.1747524569396468   max=0.8444302643579830   mean=-0.2719212961752963
+    ssh   min=-2.1845582485095099   max=0.8341590228046922   mean=-0.2740380428121212
     mld   min=2.4370615877370523   max=4593.1523423819935488   mean=193.6208071648662212
 layer_depth   min=2.4370615877370523   max=5587.8978052886013757   mean=1243.8222388394297013
 
@@ -23,10 +23,10 @@ layer_depth   min=2.4370615877370523   max=5587.8978052886013757   mean=1243.822
 
 
 
-CostJb   : Nonlinear Jb = 48.1596997487495315
-CostJo   : Nonlinear Jo(SeaSurfaceTemp) = 489.5537720365405221, nobs = 123, Jo/n = 3.9801119677767520, err = 0.3706289928395726
-CostJo   : Nonlinear Jo(SeaSurfaceSalinity) = 60.3729231990664914, nobs = 50, Jo/n = 1.2074584639813297, err = 1.0000000000000000
-CostJo   : Nonlinear Jo(ADT) = 153.7954077362375926, nobs = 91, Jo/n = 1.6900594256729407, err = 0.1000000014901161
-CostJo   : Nonlinear Jo(InsituTemperature) = 297.7986374325980705, nobs = 203, Jo/n = 1.4669883617369364, err = 0.9003481385241890
-CostJo   : Nonlinear Jo(InsituSalinity) = 146.4616360198882035, nobs = 218, Jo/n = 0.6718423670637074, err = 0.6081969954577753
-CostFunction: Nonlinear J = 1196.1420761730803406
+CostJb   : Nonlinear Jb = 52.3885394491905743
+CostJo   : Nonlinear Jo(SeaSurfaceTemp) = 504.6395987478990719, nobs = 123, Jo/n = 4.1027609654300736, err = 0.3709896629900664
+CostJo   : Nonlinear Jo(SeaSurfaceSalinity) = 61.4804594017453780, nobs = 50, Jo/n = 1.2296091880349076, err = 1.0000000000000000
+CostJo   : Nonlinear Jo(ADT) = 163.5984040081233388, nobs = 93, Jo/n = 1.7591226237432618, err = 0.1000000014901161
+CostJo   : Nonlinear Jo(InsituTemperature) = 343.9723498786432856, nobs = 205, Jo/n = 1.6779139018470404, err = 0.9016910826690303
+CostJo   : Nonlinear Jo(InsituSalinity) = 150.1449136404517617, nobs = 218, Jo/n = 0.6887381359653750, err = 0.6081969954577753
+CostFunction: Nonlinear J = 1276.2242651260532966
diff --git a/test/testref/3dvarfgat.test b/test/testref/3dvarfgat.test
index a4ce5339f..e303bcc5d 100644
--- a/test/testref/3dvarfgat.test
+++ b/test/testref/3dvarfgat.test
@@ -1,23 +1,23 @@
 CostJb   : Nonlinear Jb = 0.0000000000000000e+00
-CostJo   : Nonlinear Jo(CoolSkin) = 4.9634479411817701e+03, nobs = 47, Jo/n = 1.0560527534429298e+02, err = 3.8763169345842458e-01
-CostJo   : Nonlinear Jo(SeaSurfaceTemp) = 3.1780656884217342e+02, nobs = 43, Jo/n = 7.3908504381900793e+00, err = 3.8581534929661854e-01
-CostJo   : Nonlinear Jo(SeaSurfaceSalinity) = 2.6167366978835638e+00, nobs = 17, Jo/n = 1.5392568811079788e-01, err = 1.0000000000000000e+00
-CostJo   : Nonlinear Jo(ADT) = 8.7166288646666715e+01, nobs = 29, Jo/n = 3.0057340912643693e+00, err = 1.0000000149011612e-01
-CostJo   : Nonlinear Jo(InsituTemperature) = 3.7309178313712394e+01, nobs = 31, Jo/n = 1.2035218810874966e+00, err = 9.0893991658246809e-01
-CostJo   : Nonlinear Jo(InsituSalinity) = 1.6779590930092294e+01, nobs = 33, Jo/n = 5.0847245242703920e-01, err = 6.1042996833236707e-01
-CostJo   : Nonlinear Jo(SeaIceFraction) = 1.9943705186059643e+02, nobs = 24, Jo/n = 8.3098771608581838e+00, err = 1.0000000149011612e-01
-CostJo   : Nonlinear Jo(SurfaceU) = 4.6119926561840613e-01, nobs = 47, Jo/n = 9.8127503323065142e-03, err = 3.8763169345842458e-01
-CostFunction: Nonlinear J = 5.6250245557385133e+03
-RPCGMinimizer: reduction in residual norm = 5.5851410651771616e-01
+CostJo   : Nonlinear Jo(CoolSkin) = 4.9663191368230073e+03, nobs = 47, Jo/n = 1.0566636461325547e+02, err = 3.8763169345842458e-01
+CostJo   : Nonlinear Jo(SeaSurfaceTemp) = 3.1454830650364431e+02, nobs = 43, Jo/n = 7.3150768954335890e+00, err = 3.8581534929661854e-01
+CostJo   : Nonlinear Jo(SeaSurfaceSalinity) = 2.6796945834674828e+00, nobs = 17, Jo/n = 1.5762909314514606e-01, err = 1.0000000000000000e+00
+CostJo   : Nonlinear Jo(ADT) = 8.7431459350288023e+01, nobs = 29, Jo/n = 3.0148779086306217e+00, err = 1.0000000149011612e-01
+CostJo   : Nonlinear Jo(InsituTemperature) = 3.7253231588255915e+01, nobs = 31, Jo/n = 1.2017171480082554e+00, err = 9.0893991658246809e-01
+CostJo   : Nonlinear Jo(InsituSalinity) = 1.6517158357758575e+01, nobs = 33, Jo/n = 5.0051995023510831e-01, err = 6.1042996833236707e-01
+CostJo   : Nonlinear Jo(SeaIceFraction) = 1.9774222033107372e+02, nobs = 24, Jo/n = 8.2392591804614046e+00, err = 1.0000000149011612e-01
+CostJo   : Nonlinear Jo(SurfaceU) = 4.8960817822670860e-01, nobs = 47, Jo/n = 1.0417195281419331e-02, err = 3.8763169345842458e-01
+CostFunction: Nonlinear J = 5.6229808157157213e+03
+RPCGMinimizer: reduction in residual norm = 5.3281254594664873e-01
 CostFunction::addIncrement: Analysis: 
   Valid time: 2018-04-15T00:00:00Z
-  cicen   min=-0.0023105442480374   max=1.0069524378029397   mean=0.1182930596900565
+  cicen   min=-0.0025196085152986   max=1.0071706445543214   mean=0.1183057762690204
   hicen   min=0.0000000000000000   max=4.0326673084246947   mean=0.4712515705773916
-   socn   min=10.7210460395083924   max=40.4416591897031168   mean=34.5448184303884460
-   tocn   min=-1.8883899372702533   max=34.9320179814565819   mean=6.0312404953269336
-   uocn   min=-0.8581725098091276   max=0.7001202371455181   mean=-0.0002591135096507
+   socn   min=10.7210460395083924   max=40.4416591897031168   mean=34.5448198880093003
+   tocn   min=-1.8883899372702533   max=34.8892618778720873   mean=6.0320785119669624
+   uocn   min=-0.8581746643255235   max=0.7001206137637456   mean=-0.0002590969203429
    vocn   min=-0.7661101215480253   max=1.4377766409421606   mean=0.0021972630399530
-    ssh   min=-1.9232051438633977   max=0.9259823164020826   mean=-0.2767178501193626
+    ssh   min=-1.9227877278789181   max=0.9259686741229282   mean=-0.2766955473164547
    hocn   min=0.0009999999999977   max=1345.6400000000003274   mean=128.6280642065023017
      sw   min=0.0000000000000000   max=0.0000000000000000   mean=0.0000000000000000
     lhf   min=0.0000000000000000   max=0.0000000000000000   mean=0.0000000000000000
@@ -36,13 +36,13 @@ layer_depth   min=2.2854716757984130   max=5658.3057467114012979   mean=1200.522
 
 
 
-CostJb   : Nonlinear Jb = 290.3526002977622511
-CostJo   : Nonlinear Jo(CoolSkin) = 3431.8171901794485166, nobs = 47, Jo/n = 73.0173870250946493, err = 0.3876316934584246
-CostJo   : Nonlinear Jo(SeaSurfaceTemp) = 979.0351134119955532, nobs = 43, Jo/n = 22.7682584514417563, err = 0.3858153492966185
-CostJo   : Nonlinear Jo(SeaSurfaceSalinity) = 2.6122184426213275, nobs = 17, Jo/n = 0.1536599083894898, err = 1.0000000000000000
-CostJo   : Nonlinear Jo(ADT) = 87.2779826946484150, nobs = 29, Jo/n = 3.0095856101602902, err = 0.1000000014901161
-CostJo   : Nonlinear Jo(InsituTemperature) = 36.2139681372996591, nobs = 31, Jo/n = 1.1681925205580534, err = 0.9089399165824681
-CostJo   : Nonlinear Jo(InsituSalinity) = 16.5334752986434736, nobs = 33, Jo/n = 0.5010144029891962, err = 0.6104299683323671
-CostJo   : Nonlinear Jo(SeaIceFraction) = 190.6815481942616373, nobs = 24, Jo/n = 7.9450645080942346, err = 0.1000000014901161
-CostJo   : Nonlinear Jo(SurfaceU) = 0.4659989433909902, nobs = 47, Jo/n = 0.0099148711359785, err = 0.3876316934584246
-CostFunction: Nonlinear J = 5034.9900956000719816
+CostJb   : Nonlinear Jb = 315.6280023700833226
+CostJo   : Nonlinear Jo(CoolSkin) = 3432.3327115698070884, nobs = 47, Jo/n = 73.0283555653150387, err = 0.3876316934584246
+CostJo   : Nonlinear Jo(SeaSurfaceTemp) = 991.0677681689126075, nobs = 43, Jo/n = 23.0480876318351768, err = 0.3858153492966185
+CostJo   : Nonlinear Jo(SeaSurfaceSalinity) = 2.6750772036788657, nobs = 17, Jo/n = 0.1573574825693450, err = 1.0000000000000000
+CostJo   : Nonlinear Jo(ADT) = 87.5150111276131639, nobs = 29, Jo/n = 3.0177590044004541, err = 0.1000000014901161
+CostJo   : Nonlinear Jo(InsituTemperature) = 36.1294749185645969, nobs = 31, Jo/n = 1.1654669328569225, err = 0.9089399165824681
+CostJo   : Nonlinear Jo(InsituSalinity) = 16.2781987282261511, nobs = 33, Jo/n = 0.4932787493401864, err = 0.6104299683323671
+CostJo   : Nonlinear Jo(SeaIceFraction) = 188.8912555943718417, nobs = 24, Jo/n = 7.8704689830988270, err = 0.1000000014901161
+CostJo   : Nonlinear Jo(SurfaceU) = 0.4936754264781402, nobs = 47, Jo/n = 0.0105037324782583, err = 0.3876316934584246
+CostFunction: Nonlinear J = 5071.0111751077347435
diff --git a/test/testref/3dvarfgat_pseudo.test b/test/testref/3dvarfgat_pseudo.test
index ab0370a95..e94490b0d 100644
--- a/test/testref/3dvarfgat_pseudo.test
+++ b/test/testref/3dvarfgat_pseudo.test
@@ -1,16 +1,16 @@
 CostJb   : Nonlinear Jb = 0.0000000000000000e+00
-CostJo   : Nonlinear Jo(SeaSurfaceTemp) = 3.1780656884219644e+02, nobs = 43, Jo/n = 7.3908504381906148e+00, err = 3.8581534929661854e-01
-CostJo   : Nonlinear Jo(SeaSurfaceSalinity) = 2.6167366978836499e+00, nobs = 17, Jo/n = 1.5392568811080293e-01, err = 1.0000000000000000e+00
-CostJo   : Nonlinear Jo(ADT) = 8.7166288554315003e+01, nobs = 29, Jo/n = 3.0057340880798278e+00, err = 1.0000000149011612e-01
-CostJo   : Nonlinear Jo(InsituTemperature) = 3.7309178313712394e+01, nobs = 31, Jo/n = 1.2035218810874966e+00, err = 9.0893991658246809e-01
-CostJo   : Nonlinear Jo(InsituSalinity) = 1.6779590930092294e+01, nobs = 33, Jo/n = 5.0847245242703920e-01, err = 6.1042996833236707e-01
-CostFunction: Nonlinear J = 4.6167836333819980e+02
-RPCGMinimizer: reduction in residual norm = 1.9447522439970583e-01
+CostJo   : Nonlinear Jo(SeaSurfaceTemp) = 3.1454830650364431e+02, nobs = 43, Jo/n = 7.3150768954335890e+00, err = 3.8581534929661854e-01
+CostJo   : Nonlinear Jo(SeaSurfaceSalinity) = 2.6796945834674828e+00, nobs = 17, Jo/n = 1.5762909314514606e-01, err = 1.0000000000000000e+00
+CostJo   : Nonlinear Jo(ADT) = 8.7431459350288023e+01, nobs = 29, Jo/n = 3.0148779086306217e+00, err = 1.0000000149011612e-01
+CostJo   : Nonlinear Jo(InsituTemperature) = 3.7253231588255915e+01, nobs = 31, Jo/n = 1.2017171480082554e+00, err = 9.0893991658246809e-01
+CostJo   : Nonlinear Jo(InsituSalinity) = 1.6517158357758575e+01, nobs = 33, Jo/n = 5.0051995023510831e-01, err = 6.1042996833236707e-01
+CostFunction: Nonlinear J = 4.5842985038341430e+02
+RPCGMinimizer: reduction in residual norm = 1.7385847916791020e-01
 CostFunction::addIncrement: Analysis: 
   Valid time: 2018-04-15T00:00:00Z
-   socn   min=10.7210460395083924   max=40.4416591897031168   mean=34.5543176703830426
-   tocn   min=-1.8883899372702533   max=31.7004645720658580   mean=6.0073919144765169
-    ssh   min=-2.0744165687315084   max=1.5044961159893417   mean=-0.3523433403563111
+   socn   min=10.7210460395083924   max=40.4416591897031168   mean=34.5502087568697647
+   tocn   min=-1.8883899372702533   max=31.7004645720658580   mean=6.0080222038698352
+    ssh   min=-1.9554800289426133   max=1.5981030680627748   mean=-0.3217603882878646
    hocn   min=0.0009999999999977   max=1345.6400000000003274   mean=128.6280642065023017
    uocn   min=-0.8581693992488438   max=0.7000954286848975   mean=-0.0002591771954069
    vocn   min=-0.7661101215480253   max=1.4377766409421606   mean=0.0021972630399530
@@ -23,10 +23,10 @@ layer_depth   min=2.2854716757984130   max=5658.3057467114012979   mean=1200.522
 
 
 
-CostJb   : Nonlinear Jb = 980.6737664494047522
-CostJo   : Nonlinear Jo(SeaSurfaceTemp) = 312.0010780620916648, nobs = 43, Jo/n = 7.2558390246998066, err = 0.3858153492966185
-CostJo   : Nonlinear Jo(SeaSurfaceSalinity) = 2.6167366978836499, nobs = 17, Jo/n = 0.1539256881108029, err = 1.0000000000000000
-CostJo   : Nonlinear Jo(ADT) = 54.1946630519478489, nobs = 29, Jo/n = 1.8687814845499258, err = 0.1000000014901161
-CostJo   : Nonlinear Jo(InsituTemperature) = 37.3091783137123940, nobs = 31, Jo/n = 1.2035218810874966, err = 0.9089399165824681
-CostJo   : Nonlinear Jo(InsituSalinity) = 16.7795909300922936, nobs = 33, Jo/n = 0.5084724524270392, err = 0.6104299683323671
-CostFunction: Nonlinear J = 1403.5750135051325742
+CostJb   : Nonlinear Jb = 919.0150955744393286
+CostJo   : Nonlinear Jo(SeaSurfaceTemp) = 309.0529570232815786, nobs = 43, Jo/n = 7.1872780703088743, err = 0.3858153492966185
+CostJo   : Nonlinear Jo(SeaSurfaceSalinity) = 2.6796945834674828, nobs = 17, Jo/n = 0.1576290931451461, err = 1.0000000000000000
+CostJo   : Nonlinear Jo(ADT) = 49.3820436185993970, nobs = 29, Jo/n = 1.7028290902965308, err = 0.1000000014901161
+CostJo   : Nonlinear Jo(InsituTemperature) = 37.2532315882559146, nobs = 31, Jo/n = 1.2017171480082554, err = 0.9089399165824681
+CostJo   : Nonlinear Jo(InsituSalinity) = 16.5171583577585750, nobs = 33, Jo/n = 0.5005199502351083, err = 0.6104299683323671
+CostFunction: Nonlinear J = 1333.9001807458021176
diff --git a/test/testref/3dvarlowres_soca.test b/test/testref/3dvarlowres_soca.test
index ac0eebe82..c7d03ae94 100644
--- a/test/testref/3dvarlowres_soca.test
+++ b/test/testref/3dvarlowres_soca.test
@@ -1,16 +1,16 @@
 CostJb   : Nonlinear Jb = 0.0000000000000000e+00
-CostJo   : Nonlinear Jo(SeaSurfaceTemp) = 7.8349275156178783e+02, nobs = 148, Jo/n = 5.2938699429850526e+00, err = 3.6139743385133083e-01
-CostJo   : Nonlinear Jo(SeaSurfaceSalinity) = 4.9679661574589510e+01, nobs = 51, Jo/n = 9.7411101126646094e-01, err = 1.0000000000000000e+00
-CostJo   : Nonlinear Jo(ADT) = 1.2846395537995437e+02, nobs = 87, Jo/n = 1.4765971882753377e+00, err = 1.0000000149011612e-01
-CostJo   : Nonlinear Jo(InsituTemperature) = 3.4357440141216807e+02, nobs = 206, Jo/n = 1.6678369000590683e+00, err = 8.9914415723945484e-01
-CostJo   : Nonlinear Jo(InsituSalinity) = 3.4655673420278907e+02, nobs = 218, Jo/n = 1.5897097899210508e+00, err = 6.0819699545777528e-01
-CostFunction: Nonlinear J = 1.6517675041312889e+03
-RPCGMinimizer: reduction in residual norm = 1.0310017206421737e-01
+CostJo   : Nonlinear Jo(SeaSurfaceTemp) = 7.8439244494350282e+02, nobs = 148, Jo/n = 5.2999489523209649e+00, err = 3.6139743385133083e-01
+CostJo   : Nonlinear Jo(SeaSurfaceSalinity) = 4.7466980363687938e+01, nobs = 50, Jo/n = 9.4933960727375877e-01, err = 1.0000000000000000e+00
+CostJo   : Nonlinear Jo(ADT) = 1.3007470268433724e+02, nobs = 87, Jo/n = 1.4951115251073246e+00, err = 1.0000000149011612e-01
+CostJo   : Nonlinear Jo(InsituTemperature) = 3.5714312918688563e+02, nobs = 206, Jo/n = 1.7337045106159497e+00, err = 8.9914415723945484e-01
+CostJo   : Nonlinear Jo(InsituSalinity) = 3.4688287576640801e+02, nobs = 218, Jo/n = 1.5912058521394863e+00, err = 6.0819699545777528e-01
+CostFunction: Nonlinear J = 1.6659601329448217e+03
+RPCGMinimizer: reduction in residual norm = 9.9782677135287226e-02
 CostFunction::addIncrement: Analysis: 
   Valid time: 2018-04-15T00:00:00Z
-   socn   min=10.7212681092863438   max=40.4416939383105003   mean=34.5453518329968858
-   tocn   min=-1.9661455814068289   max=32.0208436859786119   mean=6.0265310180261373
-    ssh   min=-2.0050589780374395   max=0.9044268051331296   mean=-0.2781242905695696
+   socn   min=10.7212756423771669   max=40.4416945098767044   mean=34.5453564323479583
+   tocn   min=-2.0266830123110844   max=32.1036940001155742   mean=6.0272609036031701
+    ssh   min=-2.0038732217217592   max=0.9043564910669520   mean=-0.2779159116943226
    hocn   min=0.0009999999999977   max=1345.6400000000003274   mean=128.6280642065023017
     mld   min=2.2854716757984130   max=4593.1533423819937525   mean=192.4109073940401515
 layer_depth   min=2.2854716757984130   max=5658.3057467114012979   mean=1200.5229536158392420
@@ -21,10 +21,10 @@ layer_depth   min=2.2854716757984130   max=5658.3057467114012979   mean=1200.522
 
 
 
-CostJb   : Nonlinear Jb = 320.1890746660009199
-CostJo   : Nonlinear Jo(SeaSurfaceTemp) = 312.1680068598821549, nobs = 148, Jo/n = 2.1092432895937985, err = 0.3613974338513308
-CostJo   : Nonlinear Jo(SeaSurfaceSalinity) = 49.7201985412463401, nobs = 51, Jo/n = 0.9749058537499282, err = 1.0000000000000000
-CostJo   : Nonlinear Jo(ADT) = 91.2175770472612442, nobs = 87, Jo/n = 1.0484778970949569, err = 0.1000000014901161
-CostJo   : Nonlinear Jo(InsituTemperature) = 138.9137970137263096, nobs = 206, Jo/n = 0.6743388204549821, err = 0.8991441572394548
-CostJo   : Nonlinear Jo(InsituSalinity) = 291.4171917516000576, nobs = 218, Jo/n = 1.3367761089522938, err = 0.6081969954577753
-CostFunction: Nonlinear J = 1203.6258458797169624
+CostJb   : Nonlinear Jb = 329.0732476489210967
+CostJo   : Nonlinear Jo(SeaSurfaceTemp) = 303.8263570221579357, nobs = 148, Jo/n = 2.0528807906902564, err = 0.3613974338513308
+CostJo   : Nonlinear Jo(SeaSurfaceSalinity) = 47.5076374445173926, nobs = 50, Jo/n = 0.9501527488903478, err = 1.0000000000000000
+CostJo   : Nonlinear Jo(ADT) = 92.8506641080837767, nobs = 87, Jo/n = 1.0672490127365952, err = 0.1000000014901161
+CostJo   : Nonlinear Jo(InsituTemperature) = 139.9612611491145913, nobs = 206, Jo/n = 0.6794235978112358, err = 0.8991441572394548
+CostJo   : Nonlinear Jo(InsituSalinity) = 292.3735989580471255, nobs = 218, Jo/n = 1.3411632979726933, err = 0.6081969954577753
+CostFunction: Nonlinear J = 1205.5927663308418687
diff --git a/test/testref/4denvar.test b/test/testref/4denvar.test
index 08ae0c7c0..62755f3a9 100644
--- a/test/testref/4denvar.test
+++ b/test/testref/4denvar.test
@@ -1,11 +1,11 @@
 CostJb   : Nonlinear Jb = 0.0000000000000000e+00
-CostJo   : Nonlinear Jo(SeaSufaceTemp) = 3.4633305116094050e+02, nobs = 47, Jo/n = 7.3687883225732023e+00, err = 3.8763169345842458e-01
-CostFunction: Nonlinear J = 3.4633305116094050e+02
-DRIPCGMinimizer: reduction in residual norm = 4.7202666285926932e-02
+CostJo   : Nonlinear Jo(SeaSufaceTemp) = 3.4227236951951875e+02, nobs = 47, Jo/n = 7.2823908408408240e+00, err = 3.8763169345842458e-01
+CostFunction: Nonlinear J = 3.4227236951951875e+02
+DRIPCGMinimizer: reduction in residual norm = 2.9427640167086277e-02
 CostFunction::addIncrement: Analysis: 
   Valid time: 2018-04-15T00:00:00Z
    socn   min=10.7210460395083924   max=40.4416591897031168   mean=34.5443897261524953
-   tocn   min=-1.8883899372702533   max=31.7004645720658580   mean=6.0172879509458523
+   tocn   min=-1.8883899372702533   max=31.7004645720658580   mean=6.0173355430613267
    uocn   min=-0.8581693992488438   max=0.7000954286848975   mean=-0.0002591771954069
    vocn   min=-0.7661101215480253   max=1.4377766409421606   mean=0.0021972630399530
     ssh   min=-1.9244847628277935   max=0.9272826517867588   mean=-0.2767903423591662
@@ -13,7 +13,7 @@ CostFunction::addIncrement: Analysis:
 layer_depth   min=2.2854716757984130   max=5658.3057467114012979   mean=1200.5229536158392420
   Valid time: 2018-04-15T03:00:00Z
    socn   min=10.7210460395083924   max=40.4416591681037119   mean=34.5459713920834233
-   tocn   min=-1.8879950757380115   max=31.6932489128285653   mean=6.0153956660849053
+   tocn   min=-1.8879950757380115   max=31.6932489128285653   mean=6.0154433795887590
    uocn   min=-0.8496117748396494   max=0.6890666856088563   mean=0.0000321916028576
    vocn   min=-0.9321576238931384   max=0.9467331169326910   mean=0.0016143192702768
     ssh   min=-1.9092725357352314   max=0.9091665415750726   mean=-0.2765781634769230
@@ -21,7 +21,7 @@ layer_depth   min=2.2854716757984130   max=5658.3057467114012979   mean=1200.522
 layer_depth   min=2.2774359470406269   max=5658.2958090376268956   mean=1200.5228825300762310
   Valid time: 2018-04-15T06:00:00Z
    socn   min=10.7210460395083924   max=40.4416591598428923   mean=34.5467399790273930
-   tocn   min=-1.8874921066163295   max=31.6858695414581994   mean=6.0149051391146102
+   tocn   min=-1.8874921066163295   max=31.6858695414581994   mean=6.0149529664159935
    uocn   min=-0.8416820355059554   max=0.6786386084306094   mean=0.0002717176669510
    vocn   min=-1.1593584500375382   max=0.9487438635750088   mean=0.0003618099374769
     ssh   min=-1.9078432544437094   max=0.9075559090823697   mean=-0.2765830659083295
@@ -30,6 +30,6 @@ layer_depth   min=2.2683307266867532   max=5658.3141547773811908   mean=1200.522
 
 
 
-CostJb   : Nonlinear Jb = 62.2625395356855904
-CostJo   : Nonlinear Jo(SeaSufaceTemp) = 83.9599573855422818, nobs = 47, Jo/n = 1.7863820720328145, err = 0.3876316934584246
-CostFunction: Nonlinear J = 146.2224969212278722
+CostJb   : Nonlinear Jb = 63.2055240418158419
+CostJo   : Nonlinear Jo(SeaSufaceTemp) = 83.7000250309115188, nobs = 47, Jo/n = 1.7808515964023728, err = 0.3876316934584246
+CostFunction: Nonlinear J = 146.9055490727273536
diff --git a/test/testref/4dhybenvar.test b/test/testref/4dhybenvar.test
index 7fb8d45d8..b0205511b 100644
--- a/test/testref/4dhybenvar.test
+++ b/test/testref/4dhybenvar.test
@@ -1,11 +1,11 @@
 CostJb   : Nonlinear Jb = 0.0000000000000000e+00
-CostJo   : Nonlinear Jo(SeaSufaceTemp) = 3.4633305116094050e+02, nobs = 47, Jo/n = 7.3687883225732023e+00, err = 3.8763169345842458e-01
-CostFunction: Nonlinear J = 3.4633305116094050e+02
-DRIPCGMinimizer: reduction in residual norm = 4.6134941183791994e-02
+CostJo   : Nonlinear Jo(SeaSufaceTemp) = 3.4227236951951875e+02, nobs = 47, Jo/n = 7.2823908408408240e+00, err = 3.8763169345842458e-01
+CostFunction: Nonlinear J = 3.4227236951951875e+02
+DRIPCGMinimizer: reduction in residual norm = 3.8967484759883209e-02
 CostFunction::addIncrement: Analysis: 
   Valid time: 2018-04-15T00:00:00Z
    socn   min=10.7210460395083924   max=40.4416591897031168   mean=34.5443897261524953
-   tocn   min=-1.8883878355555379   max=31.7004645720344698   mean=6.0159257288759251
+   tocn   min=-1.8883876408043290   max=31.7004614942092680   mean=6.0159408762362077
    uocn   min=-0.8581693992488438   max=0.7000954286848975   mean=-0.0002591771954069
    vocn   min=-0.7661101215480253   max=1.4377766409421606   mean=0.0021972630399530
     ssh   min=-1.9244847628277935   max=0.9272826517867588   mean=-0.2767903423591662
@@ -13,7 +13,7 @@ CostFunction::addIncrement: Analysis:
 layer_depth   min=2.2854716757984130   max=5658.3057467114012979   mean=1200.5229536158392420
   Valid time: 2018-04-15T03:00:00Z
    socn   min=10.7210460395083924   max=40.4416591681037119   mean=34.5459713920834233
-   tocn   min=-1.8879929740232961   max=31.6932489127971770   mean=6.0140335224317027
+   tocn   min=-1.8879927792720872   max=31.6932458349719752   mean=6.0140486914996512
    uocn   min=-0.8496117748396494   max=0.6890666856088563   mean=0.0000321916028576
    vocn   min=-0.9321576238931384   max=0.9467331169326910   mean=0.0016143192702768
     ssh   min=-1.9092725357352314   max=0.9091665415750726   mean=-0.2765781634769230
@@ -21,15 +21,15 @@ layer_depth   min=2.2854716757984130   max=5658.3057467114012979   mean=1200.522
 layer_depth   min=2.2774359470406269   max=5658.2958090376268956   mean=1200.5228825300762310
   Valid time: 2018-04-15T06:00:00Z
    socn   min=10.7210460395083924   max=40.4416591598428923   mean=34.5467399790273930
-   tocn   min=-1.8874900049016143   max=31.6858695414268112   mean=6.0135431285912020
-   uocn   min=-0.8416820355059554   max=0.6786386084306091   mean=0.0002717176669510
-   vocn   min=-1.1593584500375396   max=0.9487438635750074   mean=0.0003618099374769
-    ssh   min=-1.9078432544437021   max=0.9075559090823697   mean=-0.2765830659083295
+   tocn   min=-1.8874898101504052   max=31.6858664636016094   mean=6.0135583206575083
+   uocn   min=-0.8416820355059554   max=0.6786386084306094   mean=0.0002717176669510
+   vocn   min=-1.1593584500375382   max=0.9487438635750088   mean=0.0003618099374769
+    ssh   min=-1.9078432544437094   max=0.9075559090823697   mean=-0.2765830659083295
    hocn   min=0.0008423557142235   max=1347.8173897171482167   mean=128.6280410343347285
 layer_depth   min=2.2683307266867532   max=5658.3141547773811908   mean=1200.5223352906537002
 
 
 
-CostJb   : Nonlinear Jb = 25.1558075501197536
-CostJo   : Nonlinear Jo(SeaSufaceTemp) = 14.2457107503906588, nobs = 47, Jo/n = 0.3031002287317162, err = 0.3876316934584246
-CostFunction: Nonlinear J = 39.4015183005104106
+CostJb   : Nonlinear Jb = 25.7042557933505016
+CostJo   : Nonlinear Jo(SeaSufaceTemp) = 13.9745147688376683, nobs = 47, Jo/n = 0.2973301014646312, err = 0.3876316934584246
+CostFunction: Nonlinear J = 39.6787705621881699
diff --git a/test/testref/addincrement.test b/test/testref/addincrement.test
index 0bc2d83eb..3ceb3fef6 100644
--- a/test/testref/addincrement.test
+++ b/test/testref/addincrement.test
@@ -14,15 +14,15 @@ layer_depth   min=2.2854716757984130   max=5658.3057467114012979   mean=1200.522
 Increment: 
   Valid time: 2018-04-15T00:00:00Z
   hsnon   min=0.0000000000000000   max=0.0000000000000000   mean=0.0000000000000000
-   socn   min=-0.0001044875304995   max=0.0036310038226467   mean=0.0000026969550798
-   tocn   min=-0.0915142683818291   max=0.0933406499581936   mean=0.0008145739386363
+   socn   min=-0.0001102819904892   max=0.0038751719650420   mean=0.0000028686591882
+   tocn   min=-0.0978163505349535   max=0.0746414324010661   mean=0.0008327865579380
 State plus increment: 
   Valid time: 2018-04-15T00:00:00Z
   cicen   min=0.0000000000000000   max=1.0000000000000000   mean=0.1175125761012928
   hicen   min=0.0000000000000000   max=4.0326673084246947   mean=0.4712515705773916
   hsnon   min=0.0000000000000000   max=1.2712833951042413   mean=0.0886865877527975
-   socn   min=10.7210460395083924   max=40.4416591897031168   mean=34.5443924231075883
-   tocn   min=-1.8883899372702533   max=31.7004645720658580   mean=6.0183795322919984
+   socn   min=10.7210460395083924   max=40.4416591897031168   mean=34.5443925948116970
+   tocn   min=-1.8883899372702533   max=31.7004645720658580   mean=6.0183977449112982
     ssh   min=-1.9244847628277935   max=0.9272826517867588   mean=-0.2767903423591662
    hocn   min=0.0009999999999977   max=1345.6400000000003274   mean=128.6280642065023017
    uocn   min=-0.8581693992488438   max=0.7000954286848975   mean=-0.0002591771954069
diff --git a/test/testref/checkpointmodel.test b/test/testref/checkpointmodel.test
index 73deb3942..2c24dcf5c 100644
--- a/test/testref/checkpointmodel.test
+++ b/test/testref/checkpointmodel.test
@@ -7,12 +7,12 @@ input background:
 analysis: 
 
   Valid time: 2018-04-15T00:00:00Z
-   socn   min=10.7210460395083924   max=40.4416591897031168   mean=34.5443924231075883
-   tocn   min=-1.8883899372702533   max=31.7004645720658580   mean=6.0183795322919984
+   socn   min=10.7210460395083924   max=40.4416591897031168   mean=34.5443925948116970
+   tocn   min=-1.8883899372702533   max=31.7004645720658580   mean=6.0183977449112982
    hocn   min=0.0009999999999977   max=1345.6400000000003274   mean=128.6280642065023017
 output background: 
 
   Valid time: 2018-04-15T00:00:00Z
-   socn   min=10.7210460395083924   max=38.0000000000000000   mean=34.5397728342353858
-   tocn   min=-1.8000000000000000   max=31.7004645720658580   mean=6.0185178888673949
+   socn   min=10.7210460395083924   max=38.0000000000000000   mean=34.5397730059394945
+   tocn   min=-1.8000000000000000   max=31.7004645720658580   mean=6.0185360887081467
    hocn   min=0.0009999999999977   max=1345.6400000000003274   mean=128.6280642065023017
diff --git a/test/testref/convertincrement.test b/test/testref/convertincrement.test
index beba15c44..269f7ce15 100644
--- a/test/testref/convertincrement.test
+++ b/test/testref/convertincrement.test
@@ -1,8 +1,8 @@
 Input increment: 
   Valid time: 2018-04-15T00:00:00Z
-   tocn   min=-0.0915142683818291   max=0.0933406499581936   mean=0.0008145739386363
-   socn   min=-0.0001044875304995   max=0.0036310038226467   mean=0.0000026969550798
-    ssh   min=-0.0001307321269105   max=0.0005881912327456   mean=0.0000041378484079
+   tocn   min=-0.0978163505349535   max=0.0746414324010661   mean=0.0008327865579380
+   socn   min=-0.0001102819904892   max=0.0038751719650420   mean=0.0000028686591882
+    ssh   min=-0.0001366142872018   max=0.0006271514953045   mean=0.0000040810294429
    hocn   min=0.0000000000000000   max=0.0000000000000000   mean=0.0000000000000000
 Trajectory state: 
   Valid time: 2018-04-15T00:00:00Z
@@ -14,7 +14,7 @@ layer_depth   min=2.2854716757984130   max=5658.3057467114012979   mean=1200.522
     mld   min=2.2854716757984130   max=4593.1533423819937525   mean=192.4109073940401515
 Output increment: 
   Valid time: 2018-04-15T00:00:00Z
-   tocn   min=-0.0915142683818291   max=0.0933406499581936   mean=0.0008145739386363
-   socn   min=-0.0001044875304995   max=0.0036310038226467   mean=        -nan
-    ssh   min=-0.0005470853696631   max=0.0015881623198805   mean=0.0000295468588674
+   tocn   min=-0.0978163505349535   max=0.0746414324010661   mean=0.0008327865579380
+   socn   min=-0.0001102819904892   max=0.0038751719650420   mean=        -nan
+    ssh   min=-0.0005848471049477   max=0.0016934357214616   mean=0.0000304786595822
    hocn   min=0.0000000000000000   max=0.0000000000000000   mean=0.0000000000000000
diff --git a/test/testref/convertstate_soca2cice.test b/test/testref/convertstate_soca2cice.test
index cd2ce9a05..03ce06fe1 100644
--- a/test/testref/convertstate_soca2cice.test
+++ b/test/testref/convertstate_soca2cice.test
@@ -1,9 +1,9 @@
 Input state: 
   Valid time: 2018-04-15T00:00:00Z
-   tocn   min=-1.8883899372702533   max=31.7004645720658580   mean=6.0183795324849614
-   socn   min=10.7210460395083924   max=40.4416591897031168   mean=34.5443924231105086
+   tocn   min=-1.8883899372702533   max=31.7004645720658580   mean=6.0183977449113044
+   socn   min=10.7210460395083924   max=40.4416591897031168   mean=34.5443925948116970
    hocn   min=0.0009999999999977   max=1345.6400000000003274   mean=128.6280642065018753
-  cicen   min=-0.0000705791652773   max=1.0003375266646581   mean=0.1175473153660691
+  cicen   min=-0.0000563665143372   max=1.0003307403191071   mean=0.1175499389485267
   hicen   min=0.0000000000000000   max=4.0326673084246947   mean=0.4712515705773919
   hsnon   min=0.0000000000000000   max=1.2712833951042413   mean=0.0886865877527975
 Variable transform: 
@@ -11,17 +11,17 @@ Variable change from: 6 variables: tocn, socn, hocn, cicen, hicen, hsnon
 Variable change to: 6 variables: tocn, socn, hocn, cicen, hicen, hsnon
 State after variable transform: 
   Valid time: 2018-04-15T00:00:00Z
-   tocn   min=-1.8883899372702533   max=31.7004645720658580   mean=6.0183795324849614
-   socn   min=10.7210460395083924   max=40.4416591897031168   mean=34.5443924231105086
+   tocn   min=-1.8883899372702533   max=31.7004645720658580   mean=6.0183977449113044
+   socn   min=10.7210460395083924   max=40.4416591897031168   mean=34.5443925948116970
    hocn   min=0.0009999999999977   max=1345.6400000000003274   mean=128.6280642065018753
-  cicen   min=0.0000000000000000   max=1.0000000000000002   mean=0.1168346009525868
-  hicen   min=0.0000000000000000   max=4.0326673084246947   mean=0.2834293565064365
-  hsnon   min=0.0000000000000000   max=1.2712833951042413   mean=0.0433604051377627
+  cicen   min=0.0000000000000000   max=1.0000000000000002   mean=0.1168440623576960
+  hicen   min=0.0000000000000000   max=4.0326673084246947   mean=0.2834314591631395
+  hsnon   min=0.0000000000000000   max=1.2712833951042413   mean=0.0433607663153277
 Output state: 
   Valid time: 2018-04-15T00:00:00Z
-   tocn   min=-1.8883899372702533   max=31.7004645720658580   mean=6.0183795324849614
-   socn   min=10.7210460395083924   max=40.4416591897031168   mean=34.5443924231105086
+   tocn   min=-1.8883899372702533   max=31.7004645720658580   mean=6.0183977449113044
+   socn   min=10.7210460395083924   max=40.4416591897031168   mean=34.5443925948116970
    hocn   min=0.0009999999999977   max=1345.6400000000003274   mean=128.6280642065018753
-  cicen   min=0.0000000000000000   max=1.0000000000000002   mean=0.1168346009525868
-  hicen   min=0.0000000000000000   max=4.0326673084246947   mean=0.2834293565064365
-  hsnon   min=0.0000000000000000   max=1.2712833951042413   mean=0.0433604051377627
+  cicen   min=0.0000000000000000   max=1.0000000000000002   mean=0.1168440623576960
+  hicen   min=0.0000000000000000   max=4.0326673084246947   mean=0.2834314591631395
+  hsnon   min=0.0000000000000000   max=1.2712833951042413   mean=0.0433607663153277
diff --git a/test/testref/hofx_3d.test b/test/testref/hofx_3d.test
index ccc0b1db0..8edf75240 100644
--- a/test/testref/hofx_3d.test
+++ b/test/testref/hofx_3d.test
@@ -17,28 +17,28 @@ State:
     chl   min=0.0000100000000000   max=4.6719899196170607   mean=0.1184829680403295
 dummy_atm1   min=5.0000000000000000   max=5.0000000000000000   mean=5.0000000000000000
 H(x): 
-CoolSkin nobs= 200 Min=-4.1785594844341176, Max=25.2197343303870412, RMS=20.1355481819576383
+CoolSkin nobs= 198 Min=-4.1785594844341176, Max=25.2197343303870412, RMS=20.1798065565019620
 
-SeaSurfaceTemp nobs= 200 Min=-0.7441076047036211, Max=30.6558206487590255, RMS=24.2399178481423263
+SeaSurfaceTemp nobs= 198 Min=-0.7441076047036211, Max=30.6558206487590255, RMS=24.2869807837752916
 
-SeaSurfaceSalinity nobs= 100 Min=32.8836750765081902, Max=37.0801278476553193, RMS=34.7676006134783506
+SeaSurfaceSalinity nobs= 100 Min=31.5934910820588684, Max=37.0674151811271742, RMS=34.7564821394803047
 
-ADT nobs= 99 Min=-1.2334297014633941, Max=1.4139517986919645, RMS=0.7474487596218524
+ADT nobs= 99 Min=-1.2252936028075887, Max=1.4117661236744268, RMS=0.7466272517343979
 
-InsituTS nobs= 436 Min=7.7989863451624588, Max=36.0186808894957267, RMS=28.9607174181215115
+InsituTS nobs= 436 Min=7.7989863451624588, Max=36.0186808894957267, RMS=28.9635284888697910
 
-InsituTemperature nobs= 218 Min=7.7989863451624588, Max=30.4004959665826888, RMS=21.2242834824286888
+InsituTemperature nobs= 218 Min=7.7989863451624588, Max=30.4004959665826888, RMS=21.2323110390361585
 
-InsituSalinity nobs= 218 Min=34.1391348515148039, Max=36.0186808894957267, RMS=35.0282185873344005
+InsituSalinity nobs= 218 Min=34.1391348515148039, Max=36.0186808894957267, RMS=35.0280021267666299
 
-SeaIceFraction nobs= 96 Min=0.0000000000000000, Max=0.9999999999999999, RMS=0.6969513764074570
+SeaIceFraction nobs= 95 Min=0.0000000000000000, Max=0.9999999999999999, RMS=0.6912432821158787
 
-SeaIceFreeboard nobs= 286 Min=-0.0696922177149834, Max=0.2750439383765838, RMS=0.0518468018186912
+SeaIceFreeboard nobs= 287 Min=-0.0696922177149834, Max=0.2756828688484410, RMS=0.0519997530634973
 
-Chlorophyll nobs= 137 Min=0.0029618509113541, Max=2.6146242767720782, RMS=0.3736317108789671
+Chlorophyll nobs= 135 Min=0.0029618509113541, Max=2.6146242767720782, RMS=0.3570788569205419
 
-SurfaceU nobs= 200 Min=-0.6389024077305796, Max=0.3885039618613573, RMS=0.1663702932634859
+SurfaceU nobs= 198 Min=-0.6389024077305796, Max=0.3885039618613573, RMS=0.1649990940029129
 
-SurfaceV nobs= 200 Min=-0.3381740415739801, Max=0.4480817217292362, RMS=0.1155115339081705
+SurfaceV nobs= 198 Min=-0.3381740415739801, Max=0.4480817217292362, RMS=0.1162951855034772
 
 End H(x)
diff --git a/test/testref/hofx_4d.test b/test/testref/hofx_4d.test
index 0e8b20092..3058e3daa 100644
--- a/test/testref/hofx_4d.test
+++ b/test/testref/hofx_4d.test
@@ -29,22 +29,22 @@ Final state:
      lw   min=0.0000000000000000   max=88.0360870361328125   mean=31.3955290035652546
      us   min=0.0000000000000000   max=0.0191600705253526   mean=0.0072022693699425
 H(x): 
-CoolSkin nobs= 21 Min=1.4428948751562416, Max=24.1361555774724366, RMS=19.7711636669012556
+CoolSkin nobs= 21 Min=1.5387458812807040, Max=24.1715123094076603, RMS=19.7918621451636199
 
-SeaSurfaceTemp nobs= 21 Min=5.8474277121364064, Max=29.7602622236446877, RMS=26.1545701582134633
+SeaSurfaceTemp nobs= 21 Min=5.8718670810702900, Max=29.7500148889060014, RMS=26.1595537939469196
 
-SeaSurfaceSalinity nobs= 9 Min=33.2253035608204002, Max=35.5527553210047955, RMS=34.4334800605473461
+SeaSurfaceSalinity nobs= 9 Min=33.4105735192755162, Max=35.5527553210047955, RMS=34.4634723542418868
 
-ADT nobs= 14 Min=-0.7588416066076941, Max=1.5069285902191589, RMS=1.0234497296579281
+ADT nobs= 14 Min=-0.7719489083672602, Max=1.5060812065126203, RMS=1.0251122457518520
 
 InsituTemperature nobs= 10 Min=8.0070271630434053, Max=30.4007245259000989, RMS=23.8358475396981540
 
 InsituSalinity nobs= 10 Min=34.1391616396567841, Max=35.6996397017128544, RMS=35.0158426602011090
 
-SeaIceFraction nobs= 14 Min=0.0000030348019657, Max=0.9972442090213249, RMS=0.6886532456176404
+SeaIceFraction nobs= 14 Min=0.0000030681734036, Max=0.9972442090213249, RMS=0.6907768231902665
 
-SurfaceU nobs= 21 Min=-0.5064808191553183, Max=0.0839417421563783, RMS=0.2906589228496731
+SurfaceU nobs= 21 Min=-0.4992601847182797, Max=0.0839417421563783, RMS=0.2868712638574619
 
-SurfaceV nobs= 21 Min=-0.2822470601936044, Max=0.2243200432575072, RMS=0.1317010425260674
+SurfaceV nobs= 21 Min=-0.2181718516993420, Max=0.2310073801777352, RMS=0.1233649712643844
 
 End H(x)
diff --git a/test/testref/hofx_4d_pseudo.test b/test/testref/hofx_4d_pseudo.test
index e788f2b50..f08536937 100644
--- a/test/testref/hofx_4d_pseudo.test
+++ b/test/testref/hofx_4d_pseudo.test
@@ -15,11 +15,11 @@ Final state:
    uocn   min=-0.8496117748396494   max=0.6890666856088563   mean=0.0000321916028576
    vocn   min=-0.9321576238931384   max=0.9467331169326910   mean=0.0016143192702768
 H(x): 
-SeaSurfaceTemp nobs= 21 Min=5.8474277121364064, Max=29.7602622236446877, RMS=26.1545701582134633
+SeaSurfaceTemp nobs= 21 Min=5.8718670810702900, Max=29.7500148889060014, RMS=26.1595537939469196
 
-SeaSurfaceSalinity nobs= 9 Min=33.2253035608204002, Max=35.5527553210047955, RMS=34.4334800605473461
+SeaSurfaceSalinity nobs= 9 Min=33.4105735192755162, Max=35.5527553210047955, RMS=34.4634723542418868
 
-ADT nobs= 14 Min=-0.7588416066076941, Max=1.5069285902191589, RMS=1.0234497296579281
+ADT nobs= 14 Min=-0.7719489083672602, Max=1.5060812065126203, RMS=1.0251122457518520
 
 InsituTemperature nobs= 10 Min=8.0070271630434053, Max=30.4007245259000989, RMS=23.8358475396981540
 
diff --git a/test/testref/letkf.test b/test/testref/letkf.test
index 9d173ed28..0a98d10e6 100644
--- a/test/testref/letkf.test
+++ b/test/testref/letkf.test
@@ -8,7 +8,6 @@ Initial state for member 1:
    hocn   min=0.0009999999999977   max=1345.6400000000003274   mean=128.6280642065023017
     chl   min=-0.0346456156552588   max=4.6657951350551032   mean=0.1184774650131123
    biop   min=-0.0000000010792964   max=0.0000001881576257   mean=0.0000000065411171
-
 Initial state for member 2:
   Valid time: 2018-04-15T00:00:00Z
    socn   min=10.7210460395083924   max=40.4416591897031168   mean=34.5055452687632709
@@ -19,7 +18,6 @@ Initial state for member 2:
    hocn   min=0.0009999999999977   max=1345.6400000000003274   mean=128.6280642065023017
     chl   min=-0.0289855887820164   max=4.6514396809792586   mean=0.1185532231894127
    biop   min=-0.0000000013979299   max=0.0000001885820462   mean=0.0000000065387977
-
 Initial state for member 3:
   Valid time: 2018-04-15T00:00:00Z
    socn   min=10.7210460395083924   max=40.4416591897031168   mean=34.5185040290381124
@@ -30,7 +28,6 @@ Initial state for member 3:
    hocn   min=0.0009999999999977   max=1345.6400000000003274   mean=128.6280642065023017
     chl   min=-0.0242260824668528   max=4.6851634656692536   mean=0.1184642679198311
    biop   min=-0.0000000008814691   max=0.0000001879042331   mean=0.0000000065376670
-
 Initial state for member 4:
   Valid time: 2018-04-15T00:00:00Z
    socn   min=10.7210460395083924   max=40.4416591897031168   mean=34.5332012222342186
@@ -41,41 +38,40 @@ Initial state for member 4:
    hocn   min=0.0009999999999977   max=1345.6400000000003274   mean=128.6280642065023017
     chl   min=-0.0339131334494770   max=4.6937598186554741   mean=0.1184950507291665
    biop   min=-0.0000000007072830   max=0.0000001876715923   mean=0.0000000065402892
-
 H(x) for member 1:
-SeaSurfaceTemp nobs= 200 Min=-0.7441076047036211, Max=30.2955112206964401, RMS=24.1664837132353441
+SeaSurfaceTemp nobs= 198 Min=-0.7441076047036211, Max=30.2955112206964401, RMS=24.2122987687137154
 
-ADT nobs= 99 Min=-1.1923713065292145, Max=1.3916260008199295, RMS=0.7504264800238478
+ADT nobs= 99 Min=-1.1885693259805175, Max=1.3876434460945708, RMS=0.7490713220443416
 
 
 H(x) for member 2:
-SeaSurfaceTemp nobs= 200 Min=-1.2863446879039639, Max=30.4841785881009670, RMS=24.2312506827673673
+SeaSurfaceTemp nobs= 198 Min=-1.3083143229037799, Max=30.4841785881009670, RMS=24.2705103504836721
 
-ADT nobs= 99 Min=-1.1897933830438123, Max=1.3713086760738142, RMS=0.7276861479842313
+ADT nobs= 99 Min=-1.1866978354853446, Max=1.3696369252177263, RMS=0.7273001932758628
 
 
 H(x) for member 3:
-SeaSurfaceTemp nobs= 200 Min=-0.7441076047036211, Max=30.2952580792004262, RMS=24.3474545696390692
+SeaSurfaceTemp nobs= 198 Min=-0.7441076047036211, Max=30.2633211470504264, RMS=24.3959007624468036
 
-ADT nobs= 99 Min=-1.3099837851441407, Max=1.4204194886710029, RMS=0.7461543998221180
+ADT nobs= 99 Min=-1.3056888848023624, Max=1.4173142489572348, RMS=0.7445100247416564
 
 
 H(x) for member 4:
-SeaSurfaceTemp nobs= 200 Min=-0.8409264954611300, Max=30.6045092155533602, RMS=24.1878325277770436
+SeaSurfaceTemp nobs= 198 Min=-0.8582812691598130, Max=30.6045092155533602, RMS=24.2341369155624484
 
-ADT nobs= 99 Min=-1.1117642360155990, Max=1.4200699920384396, RMS=0.7328194094163875
+ADT nobs= 99 Min=-1.1230641461548929, Max=1.4181402812731112, RMS=0.7321102696208540
 
 
 H(x) ensemble background mean: 
-SeaSurfaceTemp nobs= 200 Min=-0.7441076047036211, Max=30.3823358499043188, RMS=24.2290644633085890
+SeaSurfaceTemp nobs= 198 Min=-0.7441076047036211, Max=30.3823358499043188, RMS=24.2740131595248840
 
-ADT nobs= 99 Min=-1.2006910934514909, Max=1.4008560394007965, RMS=0.7375376522615416
+ADT nobs= 99 Min=-1.1928555077586021, Max=1.3981837253856608, RMS=0.7365479206994897
 
 
 background y - H(x): 
-SeaSurfaceTemp nobs= 200 Min=-5.7729528751167436, Max=4.0109120120212314, RMS=1.3121091098438074
+SeaSurfaceTemp nobs= 198 Min=-5.9783801761836628, Max=2.8701768824562883, RMS=1.2887996706803055
 
-ADT nobs= 99 Min=-0.5287956883999732, Max=1.0899684904517966, RMS=0.2107554631800747
+ADT nobs= 99 Min=-0.5270376309721526, Max=1.1067453550051631, RMS=0.2119704150836391
 
 
 Background mean :
@@ -88,28 +84,26 @@ Background mean :
    hocn   min=0.0009999999999977   max=1345.6400000000003274   mean=128.6280642065023017
     chl   min=-0.0135973029322532   max=4.6740395250897722   mean=0.1184975017128806
    biop   min=-0.0000000004021572   max=0.0000001872794279   mean=0.0000000065394678
-
 Analysis mean :
   Valid time: 2018-04-15T00:00:00Z
-   socn   min=10.7210460395083924   max=40.4416591897031168   mean=34.5460582533093188
-   tocn   min=-4.1420151042420219   max=39.0742940340608271   mean=6.0527804302216772
-    ssh   min=-1.9808850587334716   max=0.9356684587670925   mean=-0.2892739626094919
-   uocn   min=-0.8569747929971832   max=0.6992920028847295   mean=-0.0002609623525139
-   vocn   min=-0.7669710796528207   max=1.4343633775728784   mean=0.0022635735165409
+   socn   min=10.7210460395083924   max=40.4416591897031168   mean=34.5478912661664381
+   tocn   min=-4.1972208948987912   max=39.3290074538890124   mean=6.0512110882483627
+    ssh   min=-1.9848140030728931   max=0.9355373828745192   mean=-0.2902906736803534
+   uocn   min=-0.8569209134742887   max=0.7002508980741655   mean=-0.0002576657694545
+   vocn   min=-0.7671475190900529   max=1.4343630396853562   mean=0.0022622465474656
    hocn   min=0.0009999999999977   max=1345.6400000000003274   mean=128.6280642065023017
-    chl   min=-0.0298669042139454   max=4.6772024962675101   mean=0.1184346217818456
-   biop   min=-0.0000000009544849   max=0.0000001880204003   mean=0.0000000065358238
-
+    chl   min=-0.0300159238280435   max=4.6772866112259539   mean=0.1184279274034556
+   biop   min=-0.0000000009305103   max=0.0000001880271417   mean=0.0000000065355531
 Analysis mean increment :
   Valid time: 2018-04-15T00:00:00Z
-   socn   min=-2.5111176119167560   max=2.0210914091274503   mean=0.0201010844364589
-   tocn   min=-6.9066492864987161   max=11.0386426204930999   mean=-0.0035637080858972
-    ssh   min=-0.7525950953079572   max=0.3561891477540637   mean=-0.0169447871693975
-   uocn   min=-0.0240690528881808   max=0.0350495637061460   mean=0.0000144684996823
-   vocn   min=-0.0146415870423201   max=0.0575333116935723   mean=0.0000572225702069
+   socn   min=-2.5274225441989699   max=2.0492357937819889   mean=0.0219340972935823
+   tocn   min=-6.6118883382731770   max=11.2933560403212851   mean=-0.0051330500592141
+    ssh   min=-0.7417393217708501   max=0.3594694722350922   mean=-0.0179614982402593
+   uocn   min=-0.0227572806641118   max=0.0360576350241506   mean=0.0000177650827417
+   vocn   min=-0.0148400119857092   max=0.0589268218533487   mean=0.0000558956011316
    hocn   min=0.0000000000000000   max=0.0000000000000000   mean=0.0000000000000000
-    chl   min=-0.0858357518258490   max=0.0852125535358046   mean=-0.0000628799310351
-   biop   min=-0.0000000083973247   max=0.0000000059423443   mean=-0.0000000000036440
+    chl   min=-0.0867138970309593   max=0.0866139466486184   mean=-0.0000695743094250
+   biop   min=-0.0000000083874164   max=0.0000000060228911   mean=-0.0000000000039147
 Forecast variance :
   Valid time: 2018-04-15T00:00:00Z
    socn   min=0.0000000000000000   max=1.8239276341677120   mean=0.0465030040856202
@@ -122,49 +116,49 @@ Forecast variance :
    biop   min=0.0000000000000000   max=0.0000000000000000   mean=0.0000000000000000
 Analysis variance :
   Valid time: 2018-04-15T00:00:00Z
-   socn   min=0.0000000000000000   max=1.7162477771883140   mean=0.0414300020489708
-   tocn   min=0.0000000000000000   max=8.1408342488048682   mean=0.2329541096542200
-    ssh   min=0.0000000000000000   max=0.0544573956569277   mean=0.0029605566333293
-   uocn   min=0.0000000000000000   max=0.0001168660285223   mean=0.0000012919697748
-   vocn   min=0.0000000000000000   max=0.0001011941501286   mean=0.0000006258398684
+   socn   min=0.0000000000000000   max=1.7111809285209596   mean=0.0414269250397383
+   tocn   min=0.0000000000000000   max=8.0520942417733874   mean=0.2328887960815975
+    ssh   min=0.0000000000000000   max=0.0545480722207001   mean=0.0029611734691749
+   uocn   min=0.0000000000000000   max=0.0001169268541471   mean=0.0000012923747533
+   vocn   min=0.0000000000000000   max=0.0001011960217317   mean=0.0000006256764157
    hocn   min=0.0000000000000000   max=0.0000000000000000   mean=0.0000000000000000
-    chl   min=0.0000000000071484   max=0.0043382986032581   mean=0.0000130242093228
+    chl   min=0.0000000000071484   max=0.0043382062002196   mean=0.0000130178412296
    biop   min=0.0000000000000000   max=0.0000000000000000   mean=0.0000000000000000
 H(x) for member 1:
-SeaSurfaceTemp nobs= 200 Min=-1.2017604877912396, Max=30.4622962082671478, RMS=24.0149394365831803
+SeaSurfaceTemp nobs= 198 Min=-1.0258815612072028, Max=30.5514848768195151, RMS=24.0523339209340747
 
-ADT nobs= 99 Min=-1.2514894507247609, Max=1.4076713991590450, RMS=0.7655671941924330
+ADT nobs= 99 Min=-1.2536010323787454, Max=1.4085661430421623, RMS=0.7661637612447397
 
 
 H(x) for member 2:
-SeaSurfaceTemp nobs= 200 Min=-1.3293539570723611, Max=30.8061036309186136, RMS=24.0698618459434712
+SeaSurfaceTemp nobs= 198 Min=-1.3248841328934857, Max=30.8886835300334006, RMS=24.1011667065713020
 
-ADT nobs= 99 Min=-1.2715243724648180, Max=1.3898582112779911, RMS=0.7452627490779942
+ADT nobs= 99 Min=-1.2584089055852639, Max=1.3924977527314240, RMS=0.7459256399338591
 
 
 H(x) for member 3:
-SeaSurfaceTemp nobs= 200 Min=-1.6065211291257153, Max=30.3285104733377082, RMS=24.1704608680599122
+SeaSurfaceTemp nobs= 198 Min=-1.8162655085656818, Max=30.2746049471946712, RMS=24.2094438686005660
 
-ADT nobs= 99 Min=-1.3753653310782692, Max=1.4360654938288577, RMS=0.7622439001356500
+ADT nobs= 99 Min=-1.3688636601100521, Max=1.4374903060275335, RMS=0.7621249768668742
 
 
 H(x) for member 4:
-SeaSurfaceTemp nobs= 200 Min=-0.9165133674704848, Max=30.6634693022188465, RMS=24.0341242343249810
+SeaSurfaceTemp nobs= 198 Min=-0.9114724371520095, Max=30.5814257572882475, RMS=24.0721544530341518
 
-ADT nobs= 99 Min=-1.1906839978389554, Max=1.4376260611915415, RMS=0.7495567552237363
+ADT nobs= 99 Min=-1.1774222688943927, Max=1.4400454643882674, RMS=0.7498524817033487
 
 
 H(x) ensemble analysis mean: 
-SeaSurfaceTemp nobs= 200 Min=-0.7441076047036211, Max=30.2952436529921378, RMS=24.0692059834924414
+SeaSurfaceTemp nobs= 198 Min=-0.7441076047036211, Max=30.2471652075657644, RMS=24.1056340946481278
 
-ADT nobs= 99 Min=-1.2722657880267008, Max=1.4178052913643588, RMS=0.7542778104026348
+ADT nobs= 99 Min=-1.2622867671583469, Max=1.4196499165473466, RMS=0.7546649338950588
 
 
 analysis y - H(x): 
-SeaSurfaceTemp nobs= 200 Min=-2.6781806478141661, Max=2.8028079747317669, RMS=0.7906203680435694
+SeaSurfaceTemp nobs= 198 Min=-2.5701146752327695, Max=2.0950567050864937, RMS=0.7557496379193599
 
-ADT nobs= 99 Min=-0.5178761356349733, Max=1.0300253402999773, RMS=0.1919918062649364
+ADT nobs= 99 Min=-0.5201674962288996, Max=1.0634924104278585, RMS=0.1942965278849433
 
 
-ombg RMS: 1.0799532875562872
-oman RMS: 0.6559869323956541
+ombg RMS: 1.0593930084663352
+oman RMS: 0.6271805086200282
diff --git a/test/testref/letkf_split_observer.test b/test/testref/letkf_split_observer.test
index e1fa09b9e..7aa433be3 100644
--- a/test/testref/letkf_split_observer.test
+++ b/test/testref/letkf_split_observer.test
@@ -8,7 +8,6 @@ Initial state for member 1:
    hocn   min=0.0009999999999977   max=1345.6400000000003274   mean=128.6280642065023017
     chl   min=-0.0346456156552588   max=4.6657951350551032   mean=0.1184774650131123
    biop   min=-0.0000000010792964   max=0.0000001881576257   mean=0.0000000065411171
-
 Initial state for member 2:
   Valid time: 2018-04-15T00:00:00Z
    socn   min=10.7210460395083924   max=40.4416591897031168   mean=34.5055452687632709
@@ -19,7 +18,6 @@ Initial state for member 2:
    hocn   min=0.0009999999999977   max=1345.6400000000003274   mean=128.6280642065023017
     chl   min=-0.0289855887820164   max=4.6514396809792586   mean=0.1185532231894127
    biop   min=-0.0000000013979299   max=0.0000001885820462   mean=0.0000000065387977
-
 Initial state for member 3:
   Valid time: 2018-04-15T00:00:00Z
    socn   min=10.7210460395083924   max=40.4416591897031168   mean=34.5185040290381124
@@ -30,7 +28,6 @@ Initial state for member 3:
    hocn   min=0.0009999999999977   max=1345.6400000000003274   mean=128.6280642065023017
     chl   min=-0.0242260824668528   max=4.6851634656692536   mean=0.1184642679198311
    biop   min=-0.0000000008814691   max=0.0000001879042331   mean=0.0000000065376670
-
 Initial state for member 4:
   Valid time: 2018-04-15T00:00:00Z
    socn   min=10.7210460395083924   max=40.4416591897031168   mean=34.5332012222342186
@@ -41,40 +38,39 @@ Initial state for member 4:
    hocn   min=0.0009999999999977   max=1345.6400000000003274   mean=128.6280642065023017
     chl   min=-0.0339131334494770   max=4.6937598186554741   mean=0.1184950507291665
    biop   min=-0.0000000007072830   max=0.0000001876715923   mean=0.0000000065402892
-
 H(x) for member 1:
-SeaSurfaceTemp nobs= 200 Min=-0.7441076047036211, Max=30.2955112206964401, RMS=24.1664837132353405
+SeaSurfaceTemp nobs= 198 Min=-0.7441076047036211, Max=30.2955112206964401, RMS=24.2122987687137119
 
-ADT nobs= 99 Min=-1.1923713065292147, Max=1.3916260008199293, RMS=0.7504264800238476
+ADT nobs= 99 Min=-1.1885693259805179, Max=1.3876434460945706, RMS=0.7490713220443413
 
 
 H(x) for member 2:
-SeaSurfaceTemp nobs= 200 Min=-1.2863446879039639, Max=30.4841785881009670, RMS=24.2312506827673708
+SeaSurfaceTemp nobs= 198 Min=-1.3083143229037799, Max=30.4841785881009670, RMS=24.2705103504836721
 
-ADT nobs= 99 Min=-1.1897933830438123, Max=1.3713086760738142, RMS=0.7276861479842313
+ADT nobs= 99 Min=-1.1866978354853446, Max=1.3696369252177261, RMS=0.7273001932758628
 
 
 H(x) for member 3:
-SeaSurfaceTemp nobs= 200 Min=-0.7441076047036211, Max=30.2952580792004262, RMS=24.3474545696390692
+SeaSurfaceTemp nobs= 198 Min=-0.7441076047036211, Max=30.2633211470504264, RMS=24.3959007624468001
 
-ADT nobs= 99 Min=-1.3099837851441407, Max=1.4204194886710029, RMS=0.7461543998221180
+ADT nobs= 99 Min=-1.3056888848023624, Max=1.4173142489572350, RMS=0.7445100247416564
 
 
 H(x) for member 4:
-SeaSurfaceTemp nobs= 200 Min=-0.8409264954611300, Max=30.6045092155533602, RMS=24.1878325277770436
+SeaSurfaceTemp nobs= 198 Min=-0.8582812691598130, Max=30.6045092155533602, RMS=24.2341369155624520
 
-ADT nobs= 99 Min=-1.1117642360155993, Max=1.4200699920384396, RMS=0.7328194094163875
+ADT nobs= 99 Min=-1.1230641461548931, Max=1.4181402812731112, RMS=0.7321102696208540
 
 
 H(x) ensemble background mean: 
-SeaSurfaceTemp nobs= 200 Min=-0.7441076047036211, Max=30.3823358499043188, RMS=24.2290644633085961
+SeaSurfaceTemp nobs= 198 Min=-0.7441076047036211, Max=30.3823358499043188, RMS=24.2740131595248876
 
-ADT nobs= 99 Min=-1.2006910934514909, Max=1.4008560394007965, RMS=0.7375376522615414
+ADT nobs= 99 Min=-1.1928555077586023, Max=1.3981837253856608, RMS=0.7365479206994896
 
 
 background y - H(x): 
-SeaSurfaceTemp nobs= 200 Min=-5.7729528751167436, Max=4.0109120120212314, RMS=1.3121091098438071
+SeaSurfaceTemp nobs= 198 Min=-5.9783801761836628, Max=2.8701768824562883, RMS=1.2887996706803055
 
-ADT nobs= 99 Min=-0.5287956883999732, Max=1.0899684904517966, RMS=0.2107554631800748
+ADT nobs= 99 Min=-0.5270376309721526, Max=1.1067453550051634, RMS=0.2119704150836391
 
 
diff --git a/test/testref/letkf_split_solver.test b/test/testref/letkf_split_solver.test
index 6ef7748f8..d6fc6cb51 100644
--- a/test/testref/letkf_split_solver.test
+++ b/test/testref/letkf_split_solver.test
@@ -8,7 +8,6 @@ Initial state for member 1:
    hocn   min=0.0009999999999977   max=1345.6400000000003274   mean=128.6280642065023017
     chl   min=-0.0346456156552588   max=4.6657951350551032   mean=0.1184774650131123
    biop   min=-0.0000000010792964   max=0.0000001881576257   mean=0.0000000065411171
-
 Initial state for member 2:
   Valid time: 2018-04-15T00:00:00Z
    socn   min=10.7210460395083924   max=40.4416591897031168   mean=34.5055452687632709
@@ -19,7 +18,6 @@ Initial state for member 2:
    hocn   min=0.0009999999999977   max=1345.6400000000003274   mean=128.6280642065023017
     chl   min=-0.0289855887820164   max=4.6514396809792586   mean=0.1185532231894127
    biop   min=-0.0000000013979299   max=0.0000001885820462   mean=0.0000000065387977
-
 Initial state for member 3:
   Valid time: 2018-04-15T00:00:00Z
    socn   min=10.7210460395083924   max=40.4416591897031168   mean=34.5185040290381124
@@ -30,7 +28,6 @@ Initial state for member 3:
    hocn   min=0.0009999999999977   max=1345.6400000000003274   mean=128.6280642065023017
     chl   min=-0.0242260824668528   max=4.6851634656692536   mean=0.1184642679198311
    biop   min=-0.0000000008814691   max=0.0000001879042331   mean=0.0000000065376670
-
 Initial state for member 4:
   Valid time: 2018-04-15T00:00:00Z
    socn   min=10.7210460395083924   max=40.4416591897031168   mean=34.5332012222342186
@@ -41,41 +38,40 @@ Initial state for member 4:
    hocn   min=0.0009999999999977   max=1345.6400000000003274   mean=128.6280642065023017
     chl   min=-0.0339131334494770   max=4.6937598186554741   mean=0.1184950507291665
    biop   min=-0.0000000007072830   max=0.0000001876715923   mean=0.0000000065402892
-
 H(x) for member 1:
-SeaSurfaceTemp nobs= 200 Min=-0.7441076040267944, Max=30.2955112457275391, RMS=24.1664837135364508
+SeaSurfaceTemp nobs= 198 Min=-0.7441076040267944, Max=30.2955112457275391, RMS=24.2122987500916409
 
-ADT nobs= 99 Min=-1.1923712491989136, Max=1.3916260004043579, RMS=0.7504264789495512
+ADT nobs= 99 Min=-1.1885693073272705, Max=1.3876434564590454, RMS=0.7490713250567731
 
 
 H(x) for member 2:
-SeaSurfaceTemp nobs= 200 Min=-1.2863446474075317, Max=30.4841785430908203, RMS=24.2312506631285345
+SeaSurfaceTemp nobs= 198 Min=-1.3083143234252930, Max=30.4841785430908203, RMS=24.2705103599596796
 
-ADT nobs= 99 Min=-1.1897933483123779, Max=1.3713086843490601, RMS=0.7276861474621327
+ADT nobs= 99 Min=-1.1866978406906128, Max=1.3696368932723999, RMS=0.7273001966731646
 
 
 H(x) for member 3:
-SeaSurfaceTemp nobs= 200 Min=-0.7441076040267944, Max=30.2952575683593750, RMS=24.3474545828501725
+SeaSurfaceTemp nobs= 198 Min=-0.7441076040267944, Max=30.2633209228515625, RMS=24.3959007805595007
 
-ADT nobs= 99 Min=-1.3099837303161621, Max=1.4204194545745850, RMS=0.7461543998130861
+ADT nobs= 99 Min=-1.3056888580322266, Max=1.4173142910003662, RMS=0.7445100240993107
 
 
 H(x) for member 4:
-SeaSurfaceTemp nobs= 200 Min=-0.8409264683723450, Max=30.6045093536376953, RMS=24.1878325159965399
+SeaSurfaceTemp nobs= 198 Min=-0.8582812547683716, Max=30.6045093536376953, RMS=24.2341369346136482
 
-ADT nobs= 99 Min=-1.1117641925811768, Max=1.4200699329376221, RMS=0.7328194018518309
+ADT nobs= 99 Min=-1.1230641603469849, Max=1.4181402921676636, RMS=0.7321102695972385
 
 
 H(x) ensemble background mean: 
-SeaSurfaceTemp nobs= 200 Min=-0.7441076040267944, Max=30.3823356628417969, RMS=24.2290644588743405
+SeaSurfaceTemp nobs= 198 Min=-0.7441076040267944, Max=30.3823356628417969, RMS=24.2740131666558732
 
-ADT nobs= 99 Min=-1.2006910443305969, Max=1.4008560180664062, RMS=0.7375376499812341
+ADT nobs= 99 Min=-1.1928555071353912, Max=1.3981837332248688, RMS=0.7365479220464279
 
 
 background y - H(x): 
-SeaSurfaceTemp nobs= 200 Min=-5.7729525566101074, Max=4.0109121799468994, RMS=1.3121090797868757
+SeaSurfaceTemp nobs= 198 Min=-5.9783802032470703, Max=2.8701767921447754, RMS=1.2887996776224724
 
-ADT nobs= 99 Min=-0.5287956893444061, Max=1.0899684876203537, RMS=0.2107554632064487
+ADT nobs= 99 Min=-0.5270376354455948, Max=1.1067453622817993, RMS=0.2119704157265872
 
 
 Background mean :
@@ -88,28 +84,26 @@ Background mean :
    hocn   min=0.0009999999999977   max=1345.6400000000003274   mean=128.6280642065023017
     chl   min=-0.0135973029322532   max=4.6740395250897722   mean=0.1184975017128806
    biop   min=-0.0000000004021572   max=0.0000001872794279   mean=0.0000000065394678
-
 Analysis mean :
   Valid time: 2018-04-15T00:00:00Z
-   socn   min=10.7210460395083924   max=40.4416591897031168   mean=34.5484994850473157
-   tocn   min=-4.1420150566263079   max=39.3342720079667814   mean=6.0647518510715024
-    ssh   min=-2.0240566313569937   max=0.9356684392545753   mean=-0.2895886675031061
-   uocn   min=-0.8585951774630758   max=0.6992919999268546   mean=-0.0002578079506054
-   vocn   min=-0.7653547567785886   max=1.4343633775731606   mean=0.0022673436078005
+   socn   min=10.7210460395083924   max=40.4416591897031168   mean=34.5503984640074577
+   tocn   min=-4.1972209081537981   max=39.5677248606151792   mean=6.0633185790219688
+    ssh   min=-2.0294993501073568   max=0.9355373664810137   mean=-0.2905568340492273
+   uocn   min=-0.8588076819362294   max=0.7002508977253543   mean=-0.0002540509942861
+   vocn   min=-0.7656023849181770   max=1.4343630396850893   mean=0.0022659262939969
    hocn   min=0.0009999999999977   max=1345.6400000000003274   mean=128.6280642065023017
-    chl   min=-0.0298668983694693   max=4.6772025002347997   mean=0.1184354336172894
-   biop   min=-0.0000000009550755   max=0.0000001880204011   mean=0.0000000065360569
-
+    chl   min=-0.0300159171410080   max=4.6772866143750731   mean=0.1184290957391284
+   biop   min=-0.0000000009315800   max=0.0000001880271432   mean=0.0000000065357994
 Analysis mean increment :
   Valid time: 2018-04-15T00:00:00Z
-   socn   min=-2.5111177427918214   max=2.0408300100356769   mean=0.0225423161744710
-   tocn   min=-6.8391128160208972   max=11.2986205943990541   mean=0.0084077127639269
-    ssh   min=-0.7525951148575687   max=0.3561893016140887   mean=-0.0172594920630121
-   uocn   min=-0.0238958861412669   max=0.0351143113857929   mean=0.0000176229015907
-   vocn   min=-0.0162364960804651   max=0.0600637187967529   mean=0.0000609926614664
+   socn   min=-2.5274238274785006   max=2.0671937826357407   mean=0.0244412951346048
+   tocn   min=-6.5033672400200508   max=11.5320734470474520   mean=0.0069744407143918
+    ssh   min=-0.7417392664573934   max=0.3594697046885069   mean=-0.0182276586091335
+   uocn   min=-0.0226872528830783   max=0.0360372663357967   mean=0.0000213798579101
+   vocn   min=-0.0164197190297472   max=0.0613328320247443   mean=0.0000595753476629
    hocn   min=0.0000000000000000   max=0.0000000000000000   mean=0.0000000000000000
-    chl   min=-0.0858357451014435   max=0.0852127026306109   mean=-0.0000620680955912
-   biop   min=-0.0000000083973207   max=0.0000000059423499   mean=-0.0000000000034109
+    chl   min=-0.0867139033665234   max=0.0866140254931551   mean=-0.0000684059737523
+   biop   min=-0.0000000083874254   max=0.0000000060228961   mean=-0.0000000000036684
 Forecast variance :
   Valid time: 2018-04-15T00:00:00Z
    socn   min=0.0000000000000000   max=1.8239276341677120   mean=0.0465030040856202
@@ -122,11 +116,11 @@ Forecast variance :
    biop   min=0.0000000000000000   max=0.0000000000000000   mean=0.0000000000000000
 Analysis variance :
   Valid time: 2018-04-15T00:00:00Z
-   socn   min=0.0000000000000000   max=1.7162477268253895   mean=0.0416665311260158
-   tocn   min=0.0000000000000000   max=8.1408331330837846   mean=0.2382693678773377
-    ssh   min=0.0000000000000000   max=0.0544573955250757   mean=0.0029741311300652
-   uocn   min=0.0000000000000000   max=0.0001168663061135   mean=0.0000013103497614
-   vocn   min=0.0000000000000000   max=0.0001011941501046   mean=0.0000006321304857
+   socn   min=0.0000000000000000   max=1.7111808931083179   mean=0.0416621224808086
+   tocn   min=0.0000000000000000   max=8.0520933228739526   mean=0.2382365628859260
+    ssh   min=0.0000000000000000   max=0.0545480718986743   mean=0.0029746652530810
+   uocn   min=0.0000000000000000   max=0.0001169271324847   mean=0.0000013108109821
+   vocn   min=0.0000000000000000   max=0.0001011960216953   mean=0.0000006319773143
    hocn   min=0.0000000000000000   max=0.0000000000000000   mean=0.0000000000000000
-    chl   min=0.0000000000071484   max=0.0043382984520450   mean=0.0000130946992844
+    chl   min=0.0000000000071484   max=0.0043382060437541   mean=0.0000130880843939
    biop   min=0.0000000000000000   max=0.0000000000000000   mean=0.0000000000000000