diff --git a/src/soca/AnalyticInit/soca_analytic_mod.F90 b/src/soca/AnalyticInit/soca_analytic_mod.F90
index 1f34042d8..277394474 100644
--- a/src/soca/AnalyticInit/soca_analytic_mod.F90
+++ b/src/soca/AnalyticInit/soca_analytic_mod.F90
@@ -6,6 +6,7 @@
 !> Analytic state/geoval initialization module used for testing only
 module soca_analytic_mod
 
+use atlas_module, only : atlas_field
 use kinds, only : kind_real
 use ufo_geovals_mod, only : ufo_geovals
 use ufo_sampled_locations_mod, only : ufo_sampled_locations
@@ -57,23 +58,30 @@ subroutine soca_analytic_geovals(geovals, locs)
 !!
 !! \see soca_analytic_val
 subroutine soca_analytic_state(state)
-    class(soca_state), intent(inout) :: state
+  class(soca_state), intent(inout) :: state
 
-    integer :: i, j, f, k
-    real(kind=kind_real) :: val
+  integer :: i, j, f, k, idx
 
-    do f = 1, size(state%fields)
-        do i = state%geom%isd, state%geom%ied
-            do j = state%geom%jsd, state%geom%jed
-                do k = 1, state%fields(f)%nz
-                    val = soca_analytic_val(&
-                        state%fields(f)%name, state%fields(f)%lat(i,j), &
-                        state%fields(f)%lon(i,j), k*1.0_kind_real)
-                    state%fields(f)%val(i,j,k) = val
-                end do
-            end do
+  type(atlas_field) :: field
+  real(kind=kind_real), pointer :: fdata(:,:)
+
+  do f=1, size(state%fields)
+    field = state%afieldset%field(state%fields(f)%name)
+    call field%data(fdata)
+
+    do j = state%geom%jsc, state%geom%jec
+      do i = state%geom%isc, state%geom%iec
+        idx = state%geom%atlas_ij2idx(i,j)
+        do k = 1, field%shape(1)
+          fdata(k,idx) = soca_analytic_val(&
+            state%fields(f)%name, state%fields(f)%lat(i,j), &
+            state%fields(f)%lon(i,j), k*1.0_kind_real)
         end do
+      end do
     end do
+    call field%set_dirty()
+    call field%final()
+  end do
 end subroutine
 
 
diff --git a/src/soca/Fields/soca_fields_mod.F90 b/src/soca/Fields/soca_fields_mod.F90
index 0ed3be657..5ee233bd7 100644
--- a/src/soca/Fields/soca_fields_mod.F90
+++ b/src/soca/Fields/soca_fields_mod.F90
@@ -95,12 +95,6 @@ module soca_fields_mod
   !>\copybrief soca_field_check_congruent \see soca_field_check_congruent
   procedure :: check_congruent => soca_field_check_congruent
 
-  !>\copybrief soca_field_update_halo \see soca_field_update_halo
-  procedure :: update_halo     => soca_field_update_halo
-
-  !>\copybrief soca_field_stencil_interp \see soca_field_stencil_interp
-  procedure :: stencil_interp  => soca_field_stencil_interp
-
   !>\copybrief soca_field_fill_masked \see soca_field_fill_masked
   procedure :: fill_masked     => soca_field_fill_masked
 
@@ -130,9 +124,6 @@ module soca_fields_mod
   !> \copybrief soca_fields_create \see soca_fields_create
   procedure :: create => soca_fields_create
 
-  !> \copybrief soca_fields_copy \see soca_fields_copy
-  procedure :: copy   => soca_fields_copy
-
   !> \copybrief soca_fields_delete \see soca_fields_delete
   procedure :: delete => soca_fields_delete
 
@@ -184,9 +175,6 @@ module soca_fields_mod
   !> \name misc
   !! \{
 
-  !> \copybrief soca_fields_update_halos \see soca_fields_update_halos
-  procedure :: update_halos => soca_fields_update_halos
-
   !> \copybrief soca_fields_tohpoints \see soca_fields_tohpoints
   procedure :: tohpoints  => soca_fields_tohpoints
   !> \}
@@ -236,18 +224,6 @@ subroutine soca_field_copy(self, rhs)
 end subroutine soca_field_copy
 
 
-! ------------------------------------------------------------------------------
-!> Update the data in the halo region of the field.
-!!
-!! \relates soca_fields_mod::soca_field
-!! \todo have field keep a pointer to its relevant sections of soca_geom?
-subroutine soca_field_update_halo(self, geom)
-  class(soca_field),     intent(inout) :: self
-  type(soca_geom), pointer, intent(in) :: geom !< soca_geom from soca_fields
-
-  call mpp_update_domains(self%val, geom%Domain%mpp_domain)
-end subroutine soca_field_update_halo
-
 ! ------------------------------------------------------------------------------
 !> Make sure the two fields are the same in terms of name, size, shape.
 !!
@@ -273,8 +249,8 @@ end subroutine soca_field_check_congruent
 !!
 !! Interpolation used is inverse distance weidghted, taking into
 !! consideration the mask and using at most 6 neighbors.
-subroutine soca_field_stencil_interp(self, geom, fromto)
-  class(soca_field), intent(inout) :: self
+subroutine soca_field_stencil_interp(field, geom, fromto)
+  real(kind=kind_real), allocatable, intent(inout) :: field(:,:,:)
   class(soca_geom),    intent(in) :: geom   !< geometry
   character(len=4),     intent(in) :: fromto !< "u2h", "v2h"
 
@@ -289,7 +265,7 @@ subroutine soca_field_stencil_interp(self, geom, fromto)
   real(kind=kind_real), allocatable :: masksrc_local(:,:), maskdst_local(:,:)
 
   ! Initialize temporary arrays
-  allocate(val_tmp, mold=self%val)
+  allocate(val_tmp, mold=field)
   val_tmp = 0_kind_real
 
   ! Identify source and destination grids
@@ -319,7 +295,7 @@ subroutine soca_field_stencil_interp(self, geom, fromto)
   end select
 
   ! Interpolate
-  allocate(val(6,self%nz))
+  allocate(val(6,size(field, 3)))
   do j = geom%jsc, geom%jec
      do i = geom%isc, geom%iec
         ! destination on land, skip
@@ -334,12 +310,12 @@ subroutine soca_field_stencil_interp(self, geom, fromto)
            if (masksrc_local(ij(1,sti), ij(2,sti)) == 0_kind_real) cycle
 
            ! outcroping of layers, skip
-           if (abs(self%val(ij(1,sti), ij(2,sti),1)) > val_max) cycle
+           if (abs(field(ij(1,sti), ij(2,sti),1)) > val_max) cycle
 
            ! store the valid neighbors
            lon_src(nn) = lonsrc_local(ij(1,sti), ij(2,sti))
            lat_src(nn) = latsrc_local(ij(1,sti), ij(2,sti))
-           val(nn,:) = self%val(ij(1,sti), ij(2,sti),:)
+           val(nn,:) = field(ij(1,sti), ij(2,sti),:)
            nn = nn + 1
         end do
         nn = nn - 1
@@ -352,7 +328,7 @@ subroutine soca_field_stencil_interp(self, geom, fromto)
         end if
      end do
   end do
-  self%val = val_tmp
+  field = val_tmp
 
 end subroutine soca_field_stencil_interp
 
@@ -512,42 +488,6 @@ subroutine soca_fields_delete(self)
 end subroutine
 
 
-! ------------------------------------------------------------------------------
-!> Copy the contents of \p rhs to \p self.
-!!
-!! \p self will be initialized with the variable names in \p rhs if
-!! not already initialized.
-!!
-!! \see soca_fields_init_vars
-!! \relates soca_fields_mod::soca_fields
-subroutine soca_fields_copy(self, rhs)
-  class(soca_fields), intent(inout) :: self
-  class(soca_fields),  intent(in)    :: rhs !< fields to copy from
-
-  character(len=:), allocatable :: vars_str(:)
-  integer :: i
-  type(soca_field), pointer :: rhs_fld
-
-  ! initialize the variables based on the names in rhs
-  if (.not. allocated(self%fields)) then
-    self%geom => rhs%geom
-    allocate(character(len=1024) :: vars_str(size(rhs%fields)))
-    do i=1, size(vars_str)
-      vars_str(i) = rhs%fields(i)%name
-    end do
-    call soca_fields_init_vars(self, vars_str)
-  end if
-
-  ! copy values from rhs to self, only if the variable exists
-  !  in self
-  do i=1,size(self%fields)
-    call rhs%get(self%fields(i)%name, rhs_fld)
-    call self%fields(i)%copy(rhs_fld)
-  end do
-
-end subroutine
-
-
 ! ------------------------------------------------------------------------------
 !> Get a pointer to the soca_field with the given name.
 !!
@@ -595,32 +535,28 @@ function soca_fields_has(self, name) result(res)
 end function
 
 
-! ------------------------------------------------------------------------------
-!> Update the halo region of all fields.
-!!
-!! \relates soca_fields_mod::soca_fields
-subroutine soca_fields_update_halos(self)
-  class(soca_fields), intent(inout) :: self
-  integer :: i
-
-  do i=1,size(self%fields)
-    call self%fields(i)%update_halo(self%geom)
-  end do
-end subroutine soca_fields_update_halos
-
-
 ! ------------------------------------------------------------------------------
 !> Set the value of all fields to one.
 !!
 !! \relates soca_fields_mod::soca_fields
 subroutine soca_fields_ones(self)
   class(soca_fields), intent(inout) :: self
+  type(atlas_field) :: field
+  real(kind=kind_real), pointer :: fdata(:,:)
   integer :: i
 
+  ! TODO delete
   do i = 1, size(self%fields)
     self%fields(i)%val = 1.0_kind_real
   end do
 
+  do i = 1, self%afieldset%size()
+    field = self%afieldset%field(i)
+    call field%data(fdata)
+    fdata(:,:) = 1.0_kind_real
+    call field%final()
+  end do
+
 end subroutine soca_fields_ones
 
 
@@ -630,12 +566,22 @@ end subroutine soca_fields_ones
 !! \relates soca_fields_mod::soca_fields
 subroutine soca_fields_zeros(self)
   class(soca_fields), intent(inout) :: self
+  type(atlas_field) :: field
+  real(kind=kind_real), pointer :: fdata(:,:)
   integer :: i
 
+  ! TODO delete
   do i = 1, size(self%fields)
     self%fields(i)%val = 0.0_kind_real
   end do
 
+  do i = 1, self%afieldset%size()
+    field = self%afieldset%field(i)
+    call field%data(fdata)
+    fdata(:,:) = 0.0_kind_real
+    call field%final()
+  end do
+
 end subroutine soca_fields_zeros
 
 
@@ -1259,30 +1205,57 @@ end subroutine soca_fields_write_rst
 subroutine soca_fields_tohpoints(self)
   class(soca_fields), intent(inout) :: self !< self
 
-  integer :: i
-  real(kind=kind_real), allocatable :: val(:,:,:)
-  real(kind=kind_real), pointer :: lon_out(:,:) => null()
-  real(kind=kind_real), pointer :: lat_out(:,:) => null()
+  integer :: i,j,k,n,idx
   character(len=4) :: fromto
 
-  ! Associate lon_out and lat_out with the h-grid
-  lon_out => self%geom%lon
-  lat_out => self%geom%lat
+  type(atlas_field) :: afield
+  real(kind=kind_real), pointer :: adata(:,:)
+  real(kind=kind_real), allocatable :: fdata(:,:,:)
 
   ! Apply interpolation to all fields, when necessary
-  do i=1,size(self%fields)
+  do n=1,size(self%fields)
     ! Check if already on h-points
-    if (self%fields(i)%metadata%grid == 'h') cycle
+    if (self%fields(n)%metadata%grid == 'h') cycle
 
     ! Interpolate to different location of the stencil
-    fromto = self%fields(i)%metadata%grid//'toh'
-    call self%fields(i)%stencil_interp(self%geom, fromto)
-    call self%fields(i)%update_halo(self%geom)
+    fromto = self%fields(n)%metadata%grid//'toh'
+
+    ! convert from atlas to 3d fortran field...
+    ! because I don't want to fully refactor stencil interpolation
+    allocate(fdata(self%geom%isd:self%geom%ied, self%geom%jsd:self%geom%jed, self%fields(n)%nz))
+    afield = self%aFieldset%field(self%fields(n)%name)
+    call afield%data(adata)
+    do j=self%geom%jsc, self%geom%jec
+      do i=self%geom%isc, self%geom%iec
+        idx = self%geom%atlas_ij2idx(i,j)
+        do k=1,afield%shape(1)
+          fdata(i,j,k) = adata(k, idx)
+        end do
+      end do
+    end do
+    call mpp_update_domains(fdata, self%geom%Domain%mpp_domain)
+
+    ! interp
+    call soca_field_stencil_interp(fdata, self%geom, fromto)
+
+    !copy back to atlas
+    do j=self%geom%jsc, self%geom%jec
+      do i=self%geom%isc, self%geom%iec
+        idx = self%geom%atlas_ij2idx(i,j)
+        do k=1,afield%shape(1)
+          adata(k, idx) = fdata(i,j,k)
+        end do
+      end do
+    end do
+    deallocate(fdata)
+
+    call afield%set_dirty()
+    call afield%final()
 
     ! Update grid location to h-points
-    self%fields(i)%metadata%grid = 'h'
-    self%fields(i)%lon => self%geom%lon
-    self%fields(i)%lat => self%geom%lat
+    self%fields(n)%metadata%grid = 'h'
+    self%fields(n)%lon => self%geom%lon
+    self%fields(n)%lat => self%geom%lat
  end do
 
 end subroutine soca_fields_tohpoints
@@ -1471,11 +1444,11 @@ subroutine soca_fields_sync_from_atlas(self)
         field%val(i,j,:) = real_ptr(:, self%geom%atlas_ij2idx(i,j))
       end do
     end do
-    call field%update_halo(self%geom)
+    call mpp_update_domains(field%val, self%geom%Domain%mpp_domain)
+
   call afield%final()
   end do
 
-
   call self%update_metadata()
 end subroutine
 
diff --git a/src/soca/Geometry/soca_geom_mod.F90 b/src/soca/Geometry/soca_geom_mod.F90
index 94da3c3e6..0fc175dac 100644
--- a/src/soca/Geometry/soca_geom_mod.F90
+++ b/src/soca/Geometry/soca_geom_mod.F90
@@ -329,9 +329,9 @@ subroutine soca_geom_init_fieldset(self, f_conf, gen)
   logical,                  intent(in)  :: gen !< if true, we are doing a full init
 
   integer :: i, j, n, jz
-  type(atlas_field) :: fArea, fInterpMask, fVertCoord, fGmask, fOwned, fMaskH, fMaskU, fMaskV
+  type(atlas_field) :: fArea, fInterpMask, fVertCoord, fGmask, fOwned, fMaskH, fMaskU, fMaskV, fCos, fSin
   real(kind=kind_real), pointer :: vArea(:,:), vInterpMask(:,:), vVertCoord(:,:), &
-                                   vMaskH(:,:), vMaskU(:,:), vMaskV(:,:)
+                                   vMaskH(:,:), vMaskU(:,:), vMaskV(:,:), vCos(:,:), vSin(:,:)
   integer, pointer :: vGmask(:,:), vOwned(:,:)
 
   ! variables needed for reading in atlas fields from gridspec file
@@ -370,6 +370,14 @@ subroutine soca_geom_init_fieldset(self, f_conf, gen)
   call self%fieldset%add(fGmask)
   call fGmask%data(vGmask)
 
+  fCos = self%functionspace%create_field(name='cos_rot', kind=atlas_real(kind_real), levels=1)
+  call self%fieldset%add(fCos)
+  call fCos%data(vCos)
+
+  fSin = self%functionspace%create_field(name='sin_rot', kind=atlas_real(kind_real), levels=1)
+  call self%fieldset%add(fSin)
+  call fSin%data(vSin)
+
   ! TODO temporary, this needs to be cleaned up
   ! It's here so that answers don't change when moving the gpnorm calculations from fortran
   ! to c++
@@ -396,6 +404,8 @@ subroutine soca_geom_init_fieldset(self, f_conf, gen)
       vMaskH(1, n) = self%mask2d(i,j)
       vMaskU(1, n) = self%mask2du(i,j)
       vMaskV(1, n) = self%mask2dv(i,j)
+      vCos(1, n) = self%cos_rot(i,j)
+      vSin(1, n) = self%sin_rot(i,j)
     end do
   end do
   do jz=1,self%nzo
@@ -407,6 +417,8 @@ subroutine soca_geom_init_fieldset(self, f_conf, gen)
   call fInterpMask%halo_exchange()
   call fArea%halo_exchange()
   call fVertCoord%halo_exchange()
+  call fCos%halo_exchange()
+  call fSin%halo_exchange()
 
   ! done, cleanup
   call fArea%final()
@@ -417,6 +429,8 @@ subroutine soca_geom_init_fieldset(self, f_conf, gen)
   call fMaskH%final()
   call fMaskU%final()
   call fMaskV%final()
+  call fCos%final()
+  call fSin%final()
 
   ! -----------------------------------------------------------------------------------------------
   if (gen) then
diff --git a/src/soca/Increment/Increment.cc b/src/soca/Increment/Increment.cc
index c233ccc90..7b3c22efd 100755
--- a/src/soca/Increment/Increment.cc
+++ b/src/soca/Increment/Increment.cc
@@ -61,7 +61,7 @@ namespace soca {
 
     // same geometry, just copy and quit
     if (geom == other.geom_) {
-      soca_increment_copy_f90(toFortran(), other.toFortran());
+      util::copyFieldSet(other.fieldSet_, fieldSet_);
       return;
     }
 
@@ -100,7 +100,7 @@ namespace soca {
     : Increment(other.geom_, other.vars_, other.time_)
   {
     if (copy) {
-      soca_increment_copy_f90(toFortran(), other.toFortran());
+      util::copyFieldSet(other.fieldSet_, fieldSet_);
     }
     Log::trace() << "Increment copy-created." << std::endl;
   }
@@ -110,7 +110,7 @@ namespace soca {
   Increment::Increment(const Increment & other)
     : Increment(other.geom_, other.vars_, other.time_)
   {
-    soca_increment_copy_f90(toFortran(), other.toFortran());
+    util::copyFieldSet(other.fieldSet_, fieldSet_);
     Log::trace() << "Increment copy-created." << std::endl;
   }
 
@@ -153,7 +153,24 @@ namespace soca {
     ASSERT(geom_ == rhs.geom_);
 
     time_ = rhs.time_;
-    soca_increment_copy_f90(toFortran(), rhs.toFortran());
+
+    // copy from rhs to self, only if field exists in self
+    for (const auto & otherField : rhs.fieldSet_) {
+      if (!fieldSet_.has(otherField.name())) continue;
+
+      atlas::Field field = fieldSet_.field(otherField.name());
+      auto view = atlas::array::make_view<double, 2>(field);
+      const auto otherView = atlas::array::make_view<double, 2>(otherField);
+      for (int jnode = 0; jnode < field.shape(0); ++jnode) {
+        for (int jlevel = 0; jlevel < field.shape(1); ++jlevel) {
+          view(jnode, jlevel) = otherView(jnode, jlevel);
+        }
+      }
+
+      // If either term in the sum is out-of-date, then the result will be out-of-date
+      field.set_dirty(field.dirty() || otherField.dirty());
+    }
+
     return *this;
   }
 
diff --git a/src/soca/Increment/IncrementFortran.h b/src/soca/Increment/IncrementFortran.h
index ee10aa021..7d802c616 100644
--- a/src/soca/Increment/IncrementFortran.h
+++ b/src/soca/Increment/IncrementFortran.h
@@ -33,7 +33,6 @@ namespace soca {
                                const oops::Variables &,
                                const atlas::field::FieldSetImpl *);
     void soca_increment_delete_f90(F90flds &);
-    void soca_increment_copy_f90(const F90flds &, const F90flds &);
     void soca_increment_random_f90(const F90flds &);
     void soca_increment_dirac_f90(const F90flds &,
                               const eckit::Configuration * const &);
diff --git a/src/soca/Increment/soca_increment.interface.F90 b/src/soca/Increment/soca_increment.interface.F90
index c9d5c726e..25e041541 100644
--- a/src/soca/Increment/soca_increment.interface.F90
+++ b/src/soca/Increment/soca_increment.interface.F90
@@ -80,7 +80,6 @@ subroutine soca_increment_dirac_c(c_key_self,c_conf) bind(c,name='soca_increment
 
   call soca_increment_registry%get(c_key_self,self)
   call self%dirac(fckit_configuration(c_conf))
-  call self%sync_to_atlas()
 
 end subroutine soca_increment_dirac_c
 
@@ -94,31 +93,10 @@ subroutine soca_increment_random_c(c_key_self) bind(c,name='soca_increment_rando
 
   call soca_increment_registry%get(c_key_self,self)
   call self%random()
-  call self%sync_to_atlas()
 
 end subroutine soca_increment_random_c
 
 
-! ------------------------------------------------------------------------------
-!> C++ interface for soca_increment_mod::soca_increment version of
-!! soca_fields_mod::soca_fields::copy()
-subroutine soca_increment_copy_c(c_key_self,c_key_rhs) bind(c,name='soca_increment_copy_f90')
-  integer(c_int), intent(in) :: c_key_self
-  integer(c_int), intent(in) :: c_key_rhs
-
-  type(soca_increment), pointer :: self
-  type(soca_increment), pointer :: rhs
-
-  call soca_increment_registry%get(c_key_self,self)
-  call soca_increment_registry%get(c_key_rhs,rhs)
-  call rhs%sync_from_atlas()
-
-  call self%copy(rhs)
-  call self%sync_to_atlas()
-
-end subroutine soca_increment_copy_c
-
-
 ! ------------------------------------------------------------------------------
 !> C++ interface for soca_increment_mod::soca_increment version of
 !! soca_fields_mod::soca_fields::read()
@@ -200,7 +178,6 @@ subroutine soca_increment_horiz_scales_c(c_key_self, c_conf) &
 
   call soca_increment_registry%get(c_key_self, f_self)
   call f_self%horiz_scales(fckit_configuration(c_conf))
-  call f_self%sync_to_atlas()
 
 end subroutine soca_increment_horiz_scales_c
 
@@ -220,7 +197,6 @@ subroutine soca_increment_vert_scales_c(c_key_self, c_vert) bind(c,name='soca_in
   call soca_increment_registry%get(c_key_self, f_self)
   vert = c_vert
   call f_self%vert_scales(c_vert)
-  call f_self%sync_to_atlas()
 
 end subroutine soca_increment_vert_scales_c
 
diff --git a/src/soca/Increment/soca_increment_mod.F90 b/src/soca/Increment/soca_increment_mod.F90
index 3a0d3585d..cb5503eb3 100644
--- a/src/soca/Increment/soca_increment_mod.F90
+++ b/src/soca/Increment/soca_increment_mod.F90
@@ -71,30 +71,55 @@ subroutine soca_increment_random(self)
   integer, parameter :: rseed = 1 ! constant for reproducability of tests
     ! NOTE: random seeds are not quite working the way expected,
     !  it is only set the first time normal_distribution() is called with a seed
-  integer :: jz, i
+  integer :: i, j, k, n, idx
 
   type(soca_field), pointer :: field
+  type(atlas_field) :: afield
+  real(kind=kind_real), pointer :: fdata(:,:)
+  real(kind=kind_real), allocatable :: tmp3d(:,:,:)
+
+
+  do n = 1, self%afieldset%size()
+    afield = self%afieldset%field(n)
+    field => self%fields(n) ! TODO remove this dependency
+    call afield%data(fdata)
 
-  ! set random values
-  do i = 1, size(self%fields)
-    field => self%fields(i)
-    ! TODO remove this once increment / state are fully separated
     ! NOTE: can't randomize "hocn", testIncrementInterpAD fails
-    if (field%name == "sea_water_cell_thickness") cycle
-    call normal_distribution(field%val,  0.0_kind_real, 1.0_kind_real, rseed)
-  end do
+    if (afield%name() == "sea_water_cell_thickness") then
+      call afield%final()
+      cycle
+    end if
 
-  ! mask out land, set to zero
-  do i=1,size(self%fields)
-    field => self%fields(i)
-    if (.not. associated(field%mask) ) cycle
-    do jz=1,field%nz
-      field%val(:,:,jz) = field%val(:,:,jz) * field%mask(:,:)
+    ! allocate and fill with random values
+    ! NOTE this is done in a weird way to keep answers from changing when it was refactored
+    allocate(tmp3d(self%geom%isd:self%geom%ied, self%geom%jsd:self%geom%jed, afield%shape(1)))
+    call normal_distribution(tmp3d,  0.0_kind_real, 1.0_kind_real, rseed)
+    do j=self%geom%jsc,self%geom%jec
+      do i=self%geom%isc,self%geom%iec
+        idx = self%geom%atlas_ij2idx(i,j)
+        do k=1,afield%shape(1)
+          fdata(k,idx) = tmp3d(i,j,k)
+        end do
+      end do
     end do
+    deallocate(tmp3d)
+
+    ! mask land
+    if (associated(field%mask)) then
+      do j=self%geom%jsc,self%geom%jec
+        do i=self%geom%isc,self%geom%iec
+          idx = self%geom%atlas_ij2idx(i,j)
+          do k=1,afield%shape(1)
+            fdata(k,idx) = fdata(k,idx) * field%mask(i,j)
+          end do
+        end do
+      end do
+    end if
+
+    call afield%set_dirty()
+    call afield%final()
   end do
 
-  ! update domains
-  call self%update_halos()
 end subroutine soca_increment_random
 
 
@@ -112,7 +137,21 @@ subroutine soca_increment_dirac(self, f_conf)
   integer :: ndir,n, jz
   integer,allocatable :: ixdir(:),iydir(:),izdir(:),ifdir(:)
 
-  type(soca_field), pointer :: field
+  type(atlas_field) :: field
+  real(kind=kind_real), pointer :: fdata(:,:)
+
+  ! Define field name mapping
+  character(len=:), allocatable :: field_names(:)
+  allocate(character(len=50)::field_names(9))
+  field_names(1) = "sea_water_potential_temperature"
+  field_names(2) = "sea_water_salinity"
+  field_names(3) = "sea_surface_height_above_geoid"
+  field_names(4) = "sea_ice_area_fraction"
+  field_names(5) = "sea_ice_thickness"
+  field_names(6) = "mass_concentration_of_chlorophyll_in_sea_water"
+  field_names(7) = "molar_concentration_of_biomass_in_sea_water_in_p_units"
+  field_names(8) = "eastward_sea_water_velocity"
+  field_names(9) = "northward_sea_water_velocity"
 
   ! Get Diracs size
   ndir = f_conf%get_size("ixdir")
@@ -140,40 +179,19 @@ subroutine soca_increment_dirac(self, f_conf)
   ! Setup Diracs
   call self%zeros()
   do n=1,ndir
-      ! skip this index if not in the bounds of this PE
-      if (ixdir(n) > iec .or. ixdir(n) < isc) cycle
-      if (iydir(n) > jec .or. iydir(n) < jsc) cycle
-
-    ! TODO this list is getting long, change it so that the field name
-    ! is directly used in the yaml?
-    field => null()
-    select case(ifdir(n))
-    case (1)
-      call self%get("sea_water_potential_temperature", field)
-    case (2)
-      call self%get("sea_water_salinity", field)
-    case (3)
-      call self%get("sea_surface_height_above_geoid", field)
-    case (4)
-      call self%get("sea_ice_area_fraction", field)
-    case (5)
-      call self%get("sea_ice_thickness", field)
-    case (6)
-      call self%get("mass_concentration_of_chlorophyll_in_sea_water", field)
-    case (7)
-      call self%get("molar_concentration_of_biomass_in_sea_water_in_p_units", field)
-    case (8)
-      call self%get("eastward_sea_water_velocity", field)
-    case (9)
-      call self%get("northward_sea_water_velocity", field)
-    case default
-      ! TODO print error that out of range
-    end select
-    if (associated(field)) then
-      jz = 1
-      if (field%nz > 1) jz = izdir(n)
-      field%val(ixdir(n),iydir(n),izdir(n)) = 1.0
-    end if
+    ! skip this index if not in the bounds of this PE
+     if (ixdir(n) > iec .or. ixdir(n) < isc) cycle
+     if (iydir(n) > jec .or. iydir(n) < jsc) cycle
+
+    ! get field
+    if (ifdir(n) <= 0 .or. ifdir(n) > 9) cycle
+    field = self%afieldset%field(field_names(ifdir(n)))
+    call field%data(fdata)
+
+    ! set dirac
+    fdata(izdir(n), self%geom%atlas_ij2idx(ixdir(n),iydir(n))) = 1.0
+
+    call field%final()
   end do
 end subroutine soca_increment_dirac
 
@@ -186,26 +204,18 @@ subroutine soca_horiz_scales(self, f_conf)
   class(soca_increment),        intent(inout) :: self
   type(fckit_configuration), value, intent(in):: f_conf   !< Configuration
 
-  integer :: i, j, jz
+  integer :: n, i, j
   type(fckit_configuration) :: subconf
-  real(kind=kind_real) :: r_base, r_mult, r_min_grid, r_min, r_max
-
-  type(atlas_field) :: aField
-  real(kind=kind_real), pointer :: aFieldPtr(:,:)
-  real(kind=kind_real), allocatable :: rossby(:,:)
-
-  ! get a copy of the rossby radius from atlas
-  allocate(rossby(self%geom%isd:self%geom%ied,self%geom%jsd:self%geom%jed))
-  rossby = 0.0
-  aField = self%geom%fieldset%field("rossby_radius")
-  call aField%data(aFieldPtr)
-  do j=self%geom%jsc,self%geom%jec
-    do i=self%geom%isc,self%geom%iec
-      rossby(i,j) = aFieldPtr(1, self%geom%atlas_ij2idx(i,j))
-    end do
-  end do
-  call aField%final()
+  real(kind=kind_real) :: r_base, r_mult, r_min_grid, r_min, r_max, val
+
+  type(atlas_field) :: afield, area, rossby
+  real(kind=kind_real), pointer :: data_field(:,:), data_area(:,:), data_rossby(:,:)
 
+  ! get a copy of the input atlas fields needed
+  rossby = self%geom%fieldset%field("rossby_radius")
+  area = self%geom%fieldset%field("area")
+  call rossby%data(data_rossby)
+  call area%data(data_area)
 
   ! NOTE, this is duplicated code also present in soca_covariance_mod and possibly elsewhere.
   ! This does not belong in soca_increment_mod and should be moved out
@@ -215,29 +225,33 @@ subroutine soca_horiz_scales(self, f_conf)
   ! 2) minimum value of "min grid mult" * grid_size is imposed
   ! 3) min/max are imposed based on "min value" and "max value"
   ! 4) converted from a gaussian sigma to Gaspari-Cohn cutoff distance
-  do i=1,size(self%fields)
+  do n=1, self%afieldset%size()
+    afield = self%afieldset%field(n)
+    call afield%data(data_field)
+
     ! get parameters for correlation lengths
-    call f_conf%get_or_die(trim(self%fields(i)%name), subconf)
+    call f_conf%get_or_die(trim(afield%name()), subconf)
     if (.not. subconf%get("base value", r_base)) r_base = 0.0
     if (.not. subconf%get("rossby mult", r_mult)) r_mult = 0.0
     if (.not. subconf%get("min grid mult", r_min_grid)) r_min_grid = 1.0
     if (.not. subconf%get("min value", r_min)) r_min = 0.0
     if (.not. subconf%get("max value", r_max)) r_max = huge(r_max)
 
-    self%fields(i)%val(:,:,1) = r_base + r_mult*rossby(:,:)
-    if (r_min_grid .gt. 0.0) then
-      self%fields(i)%val(:,:,1) = max(self%fields(i)%val(:,:,1), sqrt(self%geom%cell_area)*r_min_grid)
-    end if
-    self%fields(i)%val(:,:,1) = min(r_max, self%fields(i)%val(:,:,1))
-    self%fields(i)%val(:,:,1) = max(r_min, self%fields(i)%val(:,:,1))
-    self%fields(i)%val(:,:,1) = 3.57_kind_real * self%fields(i)%val(:,:,1) ! convert from gaussian sigma to
-                                                                           ! Gaspari-Cohn half width
-
-    do jz=2,self%fields(i)%nz
-      self%fields(i)%val(:,:,jz) = self%fields(i)%val(:,:,1)
+    do i=1, afield%shape(2)
+      val = r_base + r_mult*data_rossby(1, i)
+      if (r_min_grid > 0.0) val = max(val, sqrt(data_area(1, i))*r_min_grid)
+      val = min(r_max, val)
+      val = max(r_min, val)
+      val = 3.57_kind_real * val ! convert from gaussian sigma to Gaspari-Cohn half width
+      do j=1, afield%shape(1)
+        data_field(j, i) = val
+      end do
     end do
 
+    call afield%final()
   end do
+  call rossby%final()
+  call area%final()
 end subroutine soca_horiz_scales
 
 
@@ -249,14 +263,29 @@ subroutine soca_vert_scales(self, vert)
   class(soca_increment), intent(inout) :: self
   real(kind=kind_real),  intent(in)    :: vert
 
-  integer :: i, jz
+  type(atlas_field) :: field, mask
+  real(kind=kind_real), pointer :: data_field(:,:), data_mask(:,:)
+
+  integer :: n, i, k
+
+  ! get a copy of the input atlas fields needed
+  mask = self%geom%fieldset%field("mask_h")
+  call mask%data(data_mask)
 
   ! compute scales
-  do i=1,size(self%fields)
-    do jz=1,self%fields(i)%nz
-      self%fields(i)%val(:,:,jz) = 3.57_kind_real*self%geom%mask2d(:,:)*vert
+  do n=1,self%afieldset%size()
+    field=self%afieldset%field(n)
+    call field%data(data_field)
+    do i=1,field%shape(2)
+      do k=1,field%shape(1)
+        data_field(k,i) = 3.57_kind_real * data_mask(1,i) * vert
+      end do
     end do
+
+    call field%final()
   end do
+
+  call mask%final()
 end subroutine soca_vert_scales
 ! ------------------------------------------------------------------------------
 
diff --git a/src/soca/State/State.cc b/src/soca/State/State.cc
index 6717e6cee..c039e3de7 100644
--- a/src/soca/State/State.cc
+++ b/src/soca/State/State.cc
@@ -71,7 +71,7 @@ namespace soca {
 
     // if geometry is the same, just copy and quit
     if (geom == other.geom_) {
-      soca_state_copy_f90(toFortran(), other.toFortran());
+      util::copyFieldSet(other.fieldSet_, fieldSet_);
       return;
     }
 
@@ -116,7 +116,7 @@ namespace soca {
     : Fields(other.geom_, other.vars_, other.time_)
   {
     soca_state_create_f90(keyFlds_, geom_.toFortran(), vars_, fieldSet_.get());
-    soca_state_copy_f90(toFortran(), other.toFortran());
+    util::copyFieldSet(other.fieldSet_, fieldSet_);
     Log::trace() << "State::State copied." << std::endl;
   }
 
@@ -133,7 +133,7 @@ namespace soca {
     ASSERT(geom_ == rhs.geom_);
 
     time_ = rhs.time_;
-    soca_state_copy_f90(toFortran(), rhs.toFortran());
+    util::copyFieldSet(rhs.fieldSet_, fieldSet_);
     return *this;
   }
 
@@ -142,14 +142,75 @@ namespace soca {
   // -----------------------------------------------------------------------------
   void State::rotate2north(const oops::Variables & u, const oops::Variables & v) {
     Log::trace() << "State::State rotate from logical to geographical North." << std::endl;
-    soca_state_rotate2north_f90(toFortran(), u, v);
+
+    ASSERT(u.size() == v.size());
+    for (size_t n = 0; n < u.size(); n++) {
+      const std::string & uName = u[n].name();
+      const std::string & vName = v[n].name();
+      if (!vars_.has(uName) || !vars_.has(vName)) {
+        throw eckit::UserError("State variables " + uName + " or " + vName + " not found.");
+      } else {
+        Log::info() << "rotating variables " << uName << " and " << vName << std::endl;
+      }
+
+      atlas::Field & uField = fieldSet_.field(uName);
+      atlas::Field & vField = fieldSet_.field(vName);
+      auto uView = atlas::array::make_view<double, 2>(uField);
+      auto vView = atlas::array::make_view<double, 2>(vField);
+      const auto & ghostView = atlas::array::make_view<int, 1>(uField.functionspace().ghost());
+      const auto & cosView = atlas::array::make_view<double, 2>(geom_.fields().field("cos_rot"));
+      const auto & sinView = atlas::array::make_view<double, 2>(geom_.fields().field("sin_rot"));
+
+      for (size_t i = 0; i < uField.shape(0); ++i) {
+        if (ghostView(i)) continue;
+
+        for (size_t j = 0; j < uField.shape(1); ++j) {
+          double uOrig = uView(i, j);
+          double vOrig = vView(i, j);
+          uView(i, j) = uOrig * cosView(i, 0) + vOrig * sinView(i, 0);
+          vView(i, j) = -uOrig * sinView(i, 0) + vOrig * cosView(i, 0);
+        }
+      }
+      uField.set_dirty();
+      vField.set_dirty();
+    }
   }
 
   // -----------------------------------------------------------------------------
 
   void State::rotate2grid(const oops::Variables & u, const oops::Variables & v) {
     Log::trace() << "State::State rotate from geographical to logical North." << std::endl;
-    soca_state_rotate2grid_f90(toFortran(), u, v);
+    ASSERT(u.size() == v.size());
+    for (size_t n = 0; n < u.size(); n++) {
+      const std::string & uName = u[n].name();
+      const std::string & vName = v[n].name();
+      if (!vars_.has(uName) || !vars_.has(vName)) {
+        throw eckit::UserError("State variables " + uName + " or " + vName + " not found.");
+      } else {
+        Log::info() << "rotating variables " << uName << " and " << vName << std::endl;
+      }
+
+      atlas::Field & uField = fieldSet_.field(uName);
+      atlas::Field & vField = fieldSet_.field(vName);
+      auto uView = atlas::array::make_view<double, 2>(uField);
+      auto vView = atlas::array::make_view<double, 2>(vField);
+      const auto & ghostView = atlas::array::make_view<int, 1>(uField.functionspace().ghost());
+      const auto & cosView = atlas::array::make_view<double, 2>(geom_.fields().field("cos_rot"));
+      const auto & sinView = atlas::array::make_view<double, 2>(geom_.fields().field("sin_rot"));
+
+      for (size_t i = 0; i < uField.shape(0); ++i) {
+        if (ghostView(i)) continue;
+
+        for (size_t j = 0; j < uField.shape(1); ++j) {
+          double uOrig = uView(i, j);
+          double vOrig = vView(i, j);
+          uView(i, j) = uOrig * cosView(i, 0) - vOrig * sinView(i, 0);
+          vView(i, j) = uOrig * sinView(i, 0) + vOrig * cosView(i, 0);
+        }
+      }
+      uField.set_dirty();
+      vField.set_dirty();
+    }
   }
 
   // -----------------------------------------------------------------------------
@@ -234,14 +295,48 @@ namespace soca {
 
   void State::logtrans(const oops::Variables & trvar) {
     Log::trace() << "State::State apply logarithmic transformation." << std::endl;
-    soca_state_logtrans_f90(toFortran(), trvar);
+
+    double minVal = 1.0e-6;
+    for (const auto & var : trvar) {
+      const std::string & varName = var.name();
+      if (!vars_.has(varName)) {
+        throw eckit::UserError("State variable " + varName + " not found in State.");
+      } else {
+        Log::info() << "transforming variable "  << varName << std::endl;
+      }
+
+      auto & field = fieldSet_.field(varName);
+      auto view = atlas::array::make_view<double, 2>(field);
+      for (size_t i = 0; i < field.shape(0); ++i) {
+        for (size_t j = 0; j < field.shape(1); ++j) {
+          view(i, j) = std::log(view(i, j)+minVal);
+        }
+      }
+    }
   }
 
   // -----------------------------------------------------------------------------
 
   void State::expontrans(const oops::Variables & trvar) {
     Log::trace() << "State::State apply exponential transformation." << std::endl;
-    soca_state_expontrans_f90(toFortran(), trvar);
+
+    double minVal = 1.0e-6;
+    for (const auto & var : trvar) {
+      const std::string & varName = var.name();
+      if (!vars_.has(varName)) {
+        throw eckit::UserError("State variable " + varName + " not found in State.");
+      } else {
+        Log::info() << "transforming variable "  << varName << std::endl;
+      }
+
+      auto & field = fieldSet_.field(varName);
+      auto view = atlas::array::make_view<double, 2>(field);
+      for (size_t i = 0; i < field.shape(0); ++i) {
+        for (size_t j = 0; j < field.shape(1); ++j) {
+          view(i, j) = std::exp(view(i, j))-minVal;
+        }
+      }
+    }
   }
 
   // -----------------------------------------------------------------------------
diff --git a/src/soca/State/StateFortran.h b/src/soca/State/StateFortran.h
index f80dc1eca..f1b837007 100644
--- a/src/soca/State/StateFortran.h
+++ b/src/soca/State/StateFortran.h
@@ -28,22 +28,13 @@ namespace soca {
                                const oops::Variables &,
                                const atlas::field::FieldSetImpl *);
     void soca_state_delete_f90(F90flds &);
-    void soca_state_copy_f90(const F90flds &, const F90flds &);
     void soca_state_read_file_f90(const F90flds &,
                                   const eckit::Configuration * const &,
                                   util::DateTime * const *);
     void soca_state_write_file_f90(const F90flds &,
                                    const eckit::Configuration * const &,
                                    const util::DateTime * const *);
-    void soca_state_rotate2grid_f90(const F90flds &,
-                                    const oops::Variables &,
-                                    const oops::Variables &);
-    void soca_state_rotate2north_f90(const F90flds &,
-                                     const oops::Variables &,
-                                     const oops::Variables &);
     void soca_state_tohgrid_f90(const F90flds &);
-    void soca_state_logtrans_f90(const F90flds &, const oops::Variables &);
-    void soca_state_expontrans_f90(const F90flds &, const oops::Variables &);
     void soca_state_analytic_f90(const F90flds &,
                                  const eckit::Configuration * const &,
                                  util::DateTime * const *);
diff --git a/src/soca/State/soca_state.interface.F90 b/src/soca/State/soca_state.interface.F90
index 3f78079d1..9b85f5b45 100644
--- a/src/soca/State/soca_state.interface.F90
+++ b/src/soca/State/soca_state.interface.F90
@@ -75,26 +75,6 @@ subroutine soca_state_delete_c(c_key_self) bind(c,name='soca_state_delete_f90')
 end subroutine soca_state_delete_c
 
 
-! ------------------------------------------------------------------------------
-!> C++ interface for soca_state_mod::soca_state version of
-!! soca_fields_mod::soca_fields::copy()
-subroutine soca_state_copy_c(c_key_self,c_key_rhs) bind(c,name='soca_state_copy_f90')
-    integer(c_int), intent(in) :: c_key_self
-    integer(c_int), intent(in) :: c_key_rhs
-
-    type(soca_state), pointer :: self
-    type(soca_state), pointer :: rhs
-
-    call soca_state_registry%get(c_key_self,self)
-    call soca_state_registry%get(c_key_rhs,rhs)
-
-    call rhs%sync_from_atlas()
-    call self%copy(rhs)
-    call self%sync_to_atlas()
-
-end subroutine soca_state_copy_c
-
-
 ! ------------------------------------------------------------------------------
 !> C++ interface for soca_state_mod::soca_state version of
 !! soca_fields_mod::soca_fields::read()
@@ -133,45 +113,6 @@ subroutine soca_state_write_file_c(c_key_fld, c_conf, c_dt) bind(c,name='soca_st
 end subroutine soca_state_write_file_c
 
 
-! ------------------------------------------------------------------------------
-!> C++ interface for soca_state_mod::soca_state::rotate()
-subroutine soca_state_rotate2grid_c(c_key_self, c_uvars, c_vvars) bind(c,name='soca_state_rotate2grid_f90')
-  integer(c_int), intent(in)     :: c_key_self
-  type(c_ptr), value, intent(in) :: c_uvars
-  type(c_ptr), value, intent(in) :: c_vvars
-
-  type(soca_state), pointer :: self
-  type(oops_variables) :: uvars, vvars
-
-  uvars = oops_variables(c_uvars)
-  vvars = oops_variables(c_vvars)
-
-  call soca_state_registry%get(c_key_self,self)
-  call self%rotate(coordinate="grid", uvars=uvars, vvars=vvars)
-  call self%sync_to_atlas()
-
-end subroutine soca_state_rotate2grid_c
-
-
-! ------------------------------------------------------------------------------
-!> C++ interface for soca_state_mod::soca_state::rotate()
-subroutine soca_state_rotate2north_c(c_key_self, c_uvars, c_vvars) bind(c,name='soca_state_rotate2north_f90')
-  integer(c_int),     intent(in) :: c_key_self
-  type(c_ptr), value, intent(in) :: c_uvars
-  type(c_ptr), value, intent(in) :: c_vvars
-
-  type(soca_state), pointer :: self
-  type(oops_variables) :: uvars, vvars
-
-  uvars = oops_variables(c_uvars)
-  vvars = oops_variables(c_vvars)
-
-  call soca_state_registry%get(c_key_self,self)
-  call self%rotate(coordinate="north", uvars=uvars, vvars=vvars)
-  call self%sync_to_atlas()
-
-end subroutine soca_state_rotate2north_c
-
 ! ------------------------------------------------------------------------------
 !> C++ interface for soca_state_mod::soca_state::tohgrid()
 subroutine soca_state_tohgrid_c(c_key_self) bind(c,name='soca_state_tohgrid_f90')
@@ -181,48 +122,10 @@ subroutine soca_state_tohgrid_c(c_key_self) bind(c,name='soca_state_tohgrid_f90'
 
   call soca_state_registry%get(c_key_self,self)
   call self%tohpoints()
-  call self%sync_to_atlas()
 
 end subroutine soca_state_tohgrid_c
 
 
-! ------------------------------------------------------------------------------
-!> C++ interface for soca_state_mod::soca_state::logexpon()
-subroutine soca_state_logtrans_c(c_key_self, c_trvars) bind(c,name='soca_state_logtrans_f90')
-  integer(c_int), intent(in)     :: c_key_self
-  type(c_ptr), value, intent(in) :: c_trvars
-
-  type(soca_state), pointer :: self
-  type(oops_variables) :: trvars
-
-  trvars = oops_variables(c_trvars)
-
-  call soca_state_registry%get(c_key_self,self)
-  call self%sync_from_atlas()
-  call self%logexpon(transfunc="log", trvars=trvars)
-  call self%sync_to_atlas()
-
-end subroutine soca_state_logtrans_c
-
-
-! ------------------------------------------------------------------------------
-!> C++ interface for soca_state_mod::soca_state::logexpon()
-subroutine soca_state_expontrans_c(c_key_self, c_trvars) bind(c,name='soca_state_expontrans_f90')
-  integer(c_int), intent(in)     :: c_key_self
-  type(c_ptr), value, intent(in) :: c_trvars
-
-  type(soca_state), pointer :: self
-  type(oops_variables) :: trvars
-
-  trvars = oops_variables(c_trvars)
-
-  call soca_state_registry%get(c_key_self,self)
-  call self%logexpon(transfunc="expon", trvars=trvars)
-  call self%sync_to_atlas()
-
-end subroutine soca_state_expontrans_c
-
-
 ! ------------------------------------------------------------------------------
 subroutine scoa_state_analytic_c(c_key_self, c_conf, c_dt) &
     bind(c,name='soca_state_analytic_f90')
@@ -236,7 +139,6 @@ subroutine scoa_state_analytic_c(c_key_self, c_conf, c_dt) &
   call soca_state_registry%get(c_key_self,self)
   call c_f_datetime (c_dt, fdate)
   call soca_analytic_state(self)
-  call self%sync_to_atlas()
 
 end subroutine scoa_state_analytic_c
 
diff --git a/src/soca/State/soca_state_mod.F90 b/src/soca/State/soca_state_mod.F90
index e3ab7be48..7107dc1e7 100644
--- a/src/soca/State/soca_state_mod.F90
+++ b/src/soca/State/soca_state_mod.F90
@@ -9,6 +9,7 @@ module soca_state_mod
 use logger_mod
 use kinds, only: kind_real
 use oops_variables_mod
+use atlas_module, only: atlas_field
 
 ! soca modules
 use soca_geom_mod
@@ -26,136 +27,6 @@ module soca_state_mod
 !! in the soca_fields base class
 type, public, extends(soca_fields) :: soca_state
 
-contains
-
-
-  !> \name misc
-  !! \{
-
-  !! TODO(travis) These remaning subroutines should probably be removed, and instead
-  !! live on as a non-linear variable change or saber outer block.... someday
-
-  !> \copybrief soca_state_rotate \see soca_state_rotate
-  procedure :: rotate => soca_state_rotate
-
-  !> \copybrief soca_state_logexpon \see soca_state_logexpon
-  procedure :: logexpon => soca_state_logexpon
-
-  !> \}
-
 end type
 
-
-!------------------------------------------------------------------------------
-contains
-!------------------------------------------------------------------------------
-
-
-! ------------------------------------------------------------------------------
-!> Rotate horizontal vector
-!!
-!! One or more sets of vectors, represented by corresponding u and v variables
-!! in the \p uvars and \p vvars lists are rotated to north (if \p coordinate == "north")
-!! or rotated back to the grid (if \p coordinate == "grid")
-!! \relates soca_state_mod::soca_state
-subroutine soca_state_rotate(self, coordinate, uvars, vvars)
-  class(soca_state),  intent(inout) :: self
-  character(len=*),      intent(in) :: coordinate !< "north" or "grid"
-  type(oops_variables),  intent(in) :: uvars !< list of one or more U variables
-  type(oops_variables),  intent(in) :: vvars !< list of one or more V variables
-
-  integer :: z, i
-  type(soca_field), pointer :: uocn, vocn
-  real(kind=kind_real), allocatable :: un(:,:,:), vn(:,:,:)
-  character(len=64) :: u_names, v_names
-
-  do i=1, uvars%nvars()
-    ! get (u, v) pair and make a copy
-    u_names = trim(uvars%variable(i))
-    v_names = trim(vvars%variable(i))
-    if (self%has(u_names).and.self%has(v_names)) then
-      call oops_log%info("rotating "//trim(u_names)//" "//trim(v_names))
-      call self%get(u_names, uocn)
-      call self%get(v_names, vocn)
-    else
-      ! Skip if no pair found.
-      call oops_log%info("not rotating "//trim(u_names)//" "//trim(v_names))
-      cycle
-    end if
-    allocate(un(size(uocn%val,1),size(uocn%val,2),size(uocn%val,3)))
-    allocate(vn(size(uocn%val,1),size(uocn%val,2),size(uocn%val,3)))
-    un = uocn%val
-    vn = vocn%val
-
-    select case(trim(coordinate))
-    case("north")   ! rotate (uocn, vocn) to geo north
-      do z=1,uocn%nz
-        uocn%val(:,:,z) = &
-        (self%geom%cos_rot(:,:)*un(:,:,z) + self%geom%sin_rot(:,:)*vn(:,:,z)) * uocn%mask(:,:)
-        vocn%val(:,:,z) = &
-        (- self%geom%sin_rot(:,:)*un(:,:,z) + self%geom%cos_rot(:,:)*vn(:,:,z)) * vocn%mask(:,:)
-      end do
-    case("grid")
-      do z=1,uocn%nz
-        uocn%val(:,:,z) = &
-        (self%geom%cos_rot(:,:)*un(:,:,z) - self%geom%sin_rot(:,:)*vn(:,:,z)) * uocn%mask(:,:)
-        vocn%val(:,:,z) = &
-        (self%geom%sin_rot(:,:)*un(:,:,z) + self%geom%cos_rot(:,:)*vn(:,:,z)) * vocn%mask(:,:)
-      end do
-    end select
-    deallocate(un, vn)
-
-    ! update halos
-    call uocn%update_halo(self%geom)
-    call vocn%update_halo(self%geom)
-  end do
-end subroutine soca_state_rotate
-
-
-! ------------------------------------------------------------------------------
-!> Apply logarithmic and exponential transformations
-!!
-!! \relates soca_state_mod::soca_state
-subroutine soca_state_logexpon(self, transfunc, trvars)
-  class(soca_state),  intent(inout) :: self
-  character(len=*),      intent(in) :: transfunc !< "log" or "expon"
-  type(oops_variables),  intent(in) :: trvars !< list of variables to transform
-
-  integer :: z, i
-  type(soca_field), pointer :: trocn
-  real(kind=kind_real), allocatable :: trn(:,:,:)
-  real(kind=kind_real) :: min_val = 1e-6_kind_real
-  character(len=64) :: tr_names
-
-  do i=1, trvars%nvars()
-    ! get a list variables to be transformed and make a copy
-    tr_names = trim(trvars%variable(i))
-    if (self%has(tr_names)) then
-      call oops_log%info("transforming "//trim(tr_names))
-      call self%get(tr_names, trocn)
-    else
-      ! Skip if no variable found.
-      call oops_log%info("not transforming "//trim(tr_names))
-      cycle
-    end if
-    allocate(trn(size(trocn%val,1),size(trocn%val,2),size(trocn%val,3)))
-    trn = trocn%val
-
-    select case(trim(transfunc))
-    case("log")   ! apply logarithmic transformation
-      trocn%val = log(trn + min_val)
-    case("expon") ! Apply exponential transformation
-      trocn%val = exp(trn) - min_val
-    end select
-
-    ! update halos
-    call trocn%update_halo(self%geom)
-
-    ! deallocate trn for next variable
-    deallocate(trn)
-  end do
-end subroutine soca_state_logexpon
-! ------------------------------------------------------------------------------
-
-
 end module
diff --git a/src/soca/VariableChange/Soca2Cice/Soca2Cice.F90 b/src/soca/VariableChange/Soca2Cice/Soca2Cice.F90
index 46902a6fe..bb8289770 100644
--- a/src/soca/VariableChange/Soca2Cice/Soca2Cice.F90
+++ b/src/soca/VariableChange/Soca2Cice/Soca2Cice.F90
@@ -108,9 +108,7 @@ subroutine soca_soca2cice_changevar_f90(c_key_self, c_key_geom, c_key_xin, c_key
   call soca_state_registry%get(c_key_xin, xin)
   call soca_state_registry%get(c_key_xout, xout)
 
-  call xin%sync_from_atlas()
   call self%changevar(geom, xin, xout)
-  call xout%sync_to_atlas()
 
 end subroutine
 
diff --git a/src/soca/VariableChange/Soca2Cice/soca_soca2cice_mod.F90 b/src/soca/VariableChange/Soca2Cice/soca_soca2cice_mod.F90
index c77173008..5da8ce8d7 100644
--- a/src/soca/VariableChange/Soca2Cice/soca_soca2cice_mod.F90
+++ b/src/soca/VariableChange/Soca2Cice/soca_soca2cice_mod.F90
@@ -5,7 +5,7 @@
 
 module soca_soca2cice_mod
 
-use atlas_module, only: atlas_geometry, atlas_indexkdtree
+use atlas_module, only: atlas_geometry, atlas_indexkdtree, atlas_field
 use fckit_configuration_module, only: fckit_configuration
 use fckit_exception_module, only: fckit_exception
 use fckit_mpi_module, only: fckit_mpi_comm
@@ -127,7 +127,6 @@ subroutine soca_soca2cice_changevar(self, geom, xa, xm)
 
   ! write cice restart
   call self%cice%write(geom)
-
 end subroutine soca_soca2cice_changevar
 
 ! ------------------------------------------------------------------------------
@@ -150,31 +149,37 @@ subroutine check_ice_bounds(self, geom, xm)
   type(soca_geom), target, intent(in)  :: geom
   type(soca_state),      intent(inout) :: xm
 
-  type(soca_field), pointer :: aice_ana, hice_ana, hsno_ana
+  type(atlas_field) :: aice, hice, hsno
+  real(kind=kind_real), pointer :: data_aice(:,:), data_hice(:,:), data_hsno(:,:)
 
-  ! pointers to soca fields (most likely an analysis)
-  call xm%get("sea_ice_area_fraction",aice_ana)
-  call xm%get("sea_ice_thickness",hice_ana)
-  call xm%get("sea_ice_snow_thickness",hsno_ana)
+  aice = xm%afieldset%field("sea_ice_area_fraction")
+  hice = xm%afieldset%field("sea_ice_thickness")
+  hsno = xm%afieldset%field("sea_ice_snow_thickness")
+  call aice%data(data_aice)
+  call hice%data(data_hice)
+  call hsno%data(data_hsno)
 
   ! check seaice fraction bounds
-  where (aice_ana%val<0_kind_real)
-     aice_ana%val = 0_kind_real
+  where (data_aice<0_kind_real)
+     data_aice = 0_kind_real
   end where
-  where (aice_ana%val>1_kind_real)
-     aice_ana%val = 1_kind_real
+  where (data_aice>1_kind_real)
+     data_aice = 1_kind_real
   end where
 
   ! check seaice thickness bounds
-  where (hice_ana%val<0_kind_real)
-     hice_ana%val = 0_kind_real
+  where (data_hice<0_kind_real)
+     data_hice = 0_kind_real
   end where
 
   ! check snow thickness bounds
-  where (hsno_ana%val<0_kind_real)
-     hsno_ana%val = 0_kind_real
+  where (data_hsno<0_kind_real)
+     data_hsno = 0_kind_real
   end where
 
+  call aice%final()
+  call hice%final()
+  call hsno%final()
 end subroutine check_ice_bounds
 
 ! ------------------------------------------------------------------------------
@@ -185,27 +190,30 @@ subroutine shuffle_ice(self, geom, xm)
   type(soca_geom), target, intent(in)  :: geom
   type(soca_state),      intent(inout) :: xm
 
-  real(kind=kind_real) :: aice, seaice_edge
-  integer :: i, j, k, n, ii, jj
-  type(soca_field), pointer :: s_ana, aice_ana
+  real(kind=kind_real) :: local_aice, seaice_edge
+  integer :: i, j, k, n, ii, jj, atlas_idx
   integer :: minidx(1), nn_max
   integer, allocatable :: idx(:)
   real(kind=kind_real), allocatable :: testmin(:)
   type(cice_state) :: cice_in
 
+  type(atlas_field) :: socn, aice
+  real(kind=kind_real), pointer :: data_socn(:,:), data_aice(:,:)
+
+  socn = xm%afieldset%field("sea_water_salinity")
+  aice = xm%afieldset%field("sea_ice_area_fraction")
+  call socn%data(data_socn)
+  call aice%data(data_aice)
+
   ! Make sure the search tree is smaller than the data size
   nn_max = min(self%cice%agg%n_src, self%shuffle_n)
   allocate(idx(nn_max), testmin(nn_max))
 
-  ! pointers to soca fields (most likely an analysis)
-  call xm%get("sea_water_salinity",s_ana)
-  call xm%get("sea_ice_area_fraction",aice_ana)
-
   call cice_in%copydata(self%cice)
-  do i = geom%isc, geom%iec
-     do j = geom%jsc, geom%jec
-
-        aice = aice_ana%val(i,j,1)    ! ice fraction analysis
+  do j = geom%jsc, geom%jec
+     do i = geom%isc, geom%iec
+        atlas_idx = geom%atlas_ij2idx(i,j)
+        local_aice = data_aice(1, atlas_idx)    ! ice fraction analysis
 
         ! Skip if outside of domain
         if (geom%lat(i,j)>0.0_kind_real) then
@@ -216,7 +224,7 @@ subroutine shuffle_ice(self, geom, xm)
           seaice_edge = self%antarctic%seaice_edge
         endif
         if (self%cice%aice(i,j).gt.seaice_edge) cycle     ! skip if the background has more ice than the threshold
-        if (aice.le.0.0_kind_real) then
+        if (local_aice.le.0.0_kind_real) then
            self%cice%aicen(i,j,:) = 0_kind_real
            self%cice%vicen(i,j,:) = 0_kind_real
            self%cice%vsnon(i,j,:) = 0_kind_real
@@ -226,13 +234,13 @@ subroutine shuffle_ice(self, geom, xm)
            self%cice%qice(i,j,:,:) = 0_kind_real
            self%cice%sice(i,j,:,:) = 0_kind_real
            self%cice%qsno(i,j,:,:) = 0_kind_real
-           self%cice%tsfcn(i,j,:) = liquidus_temperature_mush(s_ana%val(i,j,1))
+           self%cice%tsfcn(i,j,:) = liquidus_temperature_mush(data_socn(1, atlas_idx))
         endif
         if (self%cice%agg%n_src == 0) cycle               ! skip if there are no points on this task with ice in the background
         ! find neighbors. TODO (G): add constraint for thickness and snow depth as well
         call self%kdtree%closestPoints(geom%lon(i,j), geom%lat(i,j), nn_max, idx)
         do k = 1, nn_max
-           testmin(k) = abs(cice_in%aice(self%cice%agg%ij(1, idx(k)), self%cice%agg%ij(2, idx(k))) - aice)
+           testmin(k) = abs(cice_in%aice(self%cice%agg%ij(1, idx(k)), self%cice%agg%ij(2, idx(k))) - local_aice)
         end do
         minidx = minloc(testmin) ! I know, I rock.
         ii = self%cice%agg%ij(1, idx(minidx(1)))
@@ -258,6 +266,8 @@ subroutine shuffle_ice(self, geom, xm)
      end do
   end do
 
+  call socn%final()
+  call aice%final()
 end subroutine shuffle_ice
 
 ! ------------------------------------------------------------------------------
@@ -268,10 +278,9 @@ subroutine cleanup_ice(self, geom, xm)
   type(soca_geom), target, intent(in)  :: geom
   type(soca_state),      intent(inout) :: xm
 
-  integer :: i, j, k, ntracers
+  integer :: i, j, k, ntracers, idx
   integer :: nt_tsfc_in, nt_qice_in, nt_qsno_in, nt_sice_in
-  real(kind=kind_real) :: aice, aice0
-  type(soca_field), pointer :: t_ana, s_ana, aice_ana, hice_ana, hsno_ana
+  real(kind=kind_real) :: local_aice, aice0
   real(kind=kind_real), allocatable :: h_bounds(:)
   real(kind=kind_real), allocatable :: tracers(:,:)   ! (ntracers, ncat)
   logical, allocatable :: first_ice(:)                ! (ncat) ! For bgc and S tracers. set to true if zapping ice.
@@ -283,12 +292,20 @@ subroutine cleanup_ice(self, geom, xm)
   real(kind=kind_real), allocatable :: fiso_ocn(:)    ! isotope flux to ocean, not used as long as icepack's
                                                       ! tr_iso is false (default), so not initialized here.
 
-  ! pointers to soca fields (most likely an analysis)
-  call xm%get("sea_water_potential_temperature",t_ana)
-  call xm%get("sea_water_salinity",s_ana)
-  call xm%get("sea_ice_area_fraction",aice_ana)
-  call xm%get("sea_ice_thickness", hice_ana)
-  call xm%get("sea_ice_snow_thickness", hsno_ana)
+  type(atlas_field) :: tocn, socn, aice, hice, hsno
+  real(kind=kind_real), pointer :: data_tocn(:,:), data_socn(:,:), data_aice(:,:), data_hice(:,:), data_hsno(:,:)
+
+  ! get fields from atlas
+  tocn = xm%afieldset%field("sea_water_potential_temperature")
+  socn = xm%afieldset%field("sea_water_salinity")
+  aice = xm%afieldset%field("sea_ice_area_fraction")
+  hice = xm%afieldset%field("sea_ice_thickness")
+  hsno = xm%afieldset%field("sea_ice_snow_thickness")
+  call tocn%data(data_tocn)
+  call socn%data(data_socn)
+  call aice%data(data_aice)
+  call hice%data(data_hice)
+  call hsno%data(data_hsno)
 
   ! get thickness category bounds
   allocate(h_bounds(0:self%ncat))
@@ -337,8 +354,10 @@ subroutine cleanup_ice(self, geom, xm)
   allocate(first_ice(self%ncat))
   first_ice(:) = .true.
 
-  do i = geom%isc, geom%iec
-     do j = geom%jsc, geom%jec
+  do j = geom%jsc, geom%jec
+     do i = geom%isc, geom%iec
+        idx = geom%atlas_ij2idx(i,j)
+
         ! setup tracers at this gridpoint
         tracers(nt_tsfc,:) = self%cice%tsfcn(i,j,:)
         do k = 1, self%ice_lev
@@ -354,7 +373,7 @@ subroutine cleanup_ice(self, geom, xm)
                          h_bounds, self%cice%aicen(i,j,:), tracers, &
                          self%cice%vicen(i,j,:), self%cice%vsnon(i,j,:), &
                          ! ice and total water concentration are computed in the call using aicen
-                         aice, aice0, &
+                         local_aice, aice0, &
                          ! number of aerosol tracers, bio tracers, bio layers
                          0, 0, 0, &
                          ! aerosol flag, topo pond flag, heat capacity flag, flag for zapping ice for bgc and s tracers
@@ -376,14 +395,24 @@ subroutine cleanup_ice(self, geom, xm)
            call abor1_ftn("Soca2Cice: icepack aborted during cleanup_itd")
         endif
         ! re-compute aggregates = analysis that is effectively inserted in the restart
-        aice_ana%val(i,j,1) = sum(self%cice%aicen(i,j,:))
-        hice_ana%val(i,j,1) = sum(self%cice%vicen(i,j,:))
-        hsno_ana%val(i,j,1) = sum(self%cice%vsnon(i,j,:))
+        data_aice(1, idx) = sum(self%cice%aicen(i,j,:))
+        data_hice(1, idx) = sum(self%cice%vicen(i,j,:))
+        data_hsno(1, idx) = sum(self%cice%vsnon(i,j,:))
      end do
   end do
 
+  ! indicate dirty halos for updated fields
+  call aice%set_dirty()
+  call hice%set_dirty()
+  call hsno%set_dirty()
+
   deallocate(h_bounds, tracers, trcr_depend, trcr_base, n_trcr_strata, nt_strata, first_ice)
 
+  call tocn%final()
+  call socn%final()
+  call aice%final()
+  call hice%final()
+  call hsno%final()
 end subroutine cleanup_ice
 
 ! ------------------------------------------------------------------------------
@@ -393,17 +422,25 @@ subroutine prior_dist_rescale(self, geom, xm)
   type(soca_geom), target, intent(in)  :: geom
   type(soca_state),      intent(inout) :: xm
 
-  real(kind=kind_real) :: alpha, hice, hsno, seaice_edge, rescale_min_hice, rescale_min_hsno
-  type(soca_field), pointer :: s_ana, aice_ana, hice_ana, hsno_ana
-  integer :: c, i, j
+  real(kind=kind_real) :: alpha, local_hice, local_hsno, seaice_edge, rescale_min_hice, rescale_min_hsno
+  integer :: c, i, j, idx
+
+  type(atlas_field) :: aice, hice, hsno, socn
+  real(kind=kind_real), pointer :: data_aice(:,:), data_hice(:,:), data_hsno(:,:), data_socn(:,:)
 
-  call xm%get("sea_ice_area_fraction", aice_ana)
-  call xm%get("sea_ice_thickness", hice_ana)
-  call xm%get("sea_ice_snow_thickness", hsno_ana)
-  call xm%get("sea_water_salinity", s_ana)
+  ! get fields from atlas
+  aice = xm%afieldset%field("sea_ice_area_fraction")
+  hice = xm%afieldset%field("sea_ice_thickness")
+  hsno = xm%afieldset%field("sea_ice_snow_thickness")
+  socn = xm%afieldset%field("sea_water_salinity")
+  call aice%data(data_aice)
+  call hice%data(data_hice)
+  call hsno%data(data_hsno)
+  call socn%data(data_socn)
 
-  do i = geom%isc, geom%iec
-     do j = geom%jsc, geom%jec
+  do j = geom%jsc, geom%jec
+    do i = geom%isc, geom%iec
+      idx = geom%atlas_ij2idx(i,j)
 
         if (geom%lat(i,j)>0.0_kind_real) then
           if (.not. self%arctic%rescale_prior) cycle
@@ -419,7 +456,7 @@ subroutine prior_dist_rescale(self, geom, xm)
         if (self%cice%aice(i,j).le.seaice_edge) cycle ! Only rescale within the icepack
 
         ! rescale background to match aggregate ice concentration analysis
-        alpha = aice_ana%val(i,j,1)/self%cice%aice(i,j)
+        alpha = data_aice(1, idx)/self%cice%aice(i,j)
         self%cice%aice(i,j) = alpha * self%cice%aice(i,j)
         do c = 1, self%ncat
            self%cice%aicen(i,j,c) = alpha*self%cice%aicen(i,j,c)
@@ -428,22 +465,25 @@ subroutine prior_dist_rescale(self, geom, xm)
         end do
 
         ! adjust ice volume to match mean cell thickness
-        hice = sum(self%cice%vicen(i,j,:))
-        if (hice.gt.rescale_min_hice) then
-           alpha = hice_ana%val(i,j,1)/hice
+        local_hice = sum(self%cice%vicen(i,j,:))
+        if (local_hice.gt.rescale_min_hice) then
+           alpha = data_hice(1, idx)/local_hice
            self%cice%vicen(i,j,:) = alpha*self%cice%vicen(i,j,:)
         end if
 
         ! adjust snow volume to match mean cell thickness
-        hsno = sum(self%cice%vsnon(i,j,:))
-        if (hsno.gt.rescale_min_hsno) then
-           alpha = hsno_ana%val(i,j,1)/hsno
+        local_hsno = sum(self%cice%vsnon(i,j,:))
+        if (local_hsno.gt.rescale_min_hsno) then
+           alpha = data_hsno(1, idx)/local_hsno
            self%cice%vsnon(i,j,:) = alpha*self%cice%vsnon(i,j,:)
         end if
-
-     end do
+    end do
   end do
 
+  call aice%final()
+  call hice%final()
+  call hsno%final()
+  call socn%final()
 end subroutine prior_dist_rescale
 
 ! ------------------------------------------------------------------------------