Use the generic oops interpolator

src/soca/CMakeLists.txt

@@ -54,7 +54,6 @@ add_subdirectory(Fields)
 add_subdirectory(Geometry)
 add_subdirectory(GeometryIterator)
 add_subdirectory(Increment)
-add_subdirectory(Interpolator)
 add_subdirectory(LinearVariableChange)
 add_subdirectory(Model)
 add_subdirectory(ModelBias)

src/soca/Covariance/soca_covariance_mod.F90

@@ -280,7 +280,7 @@ subroutine soca_bump_correlation(self, horiz_convol, geom, f_conf_bump, f_conf_d
   ! 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,.true.)
+  area = pack(geom%cell_area,geom%valid_halo_mask)
   real_ptr(1,:) = area
   call afieldset%add(afield)
   call afield%final()
@@ -295,7 +295,7 @@ subroutine soca_bump_correlation(self, horiz_convol, geom, f_conf_bump, f_conf_d
   ! 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,.true.))
+  int_ptr(1,:) = int(pack(geom%mask2d,geom%valid_halo_mask))
   call afieldset%add(afield)
   call afield%final()
 
@@ -304,7 +304,7 @@ subroutine soca_bump_correlation(self, horiz_convol, geom, f_conf_bump, f_conf_d
   hmask = 0
   hmask(geom%isc:geom%iec, geom%jsc:geom%jec) = 1
   call afield%data(int_ptr)
-  int_ptr(1,:) = pack(hmask, .true.)
+  int_ptr(1,:) = pack(hmask, geom%valid_halo_mask)
   call afieldset%add(afield)
   call afield%final()
 
@@ -331,7 +331,7 @@ subroutine soca_bump_correlation(self, horiz_convol, geom, f_conf_bump, f_conf_d
     afield = geom%functionspaceInchalo%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, .true.)
+    real_ptr(1,:) = r_base + r_mult*pack(geom%rossby_radius, geom%valid_halo_mask)
     ! 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 )

src/soca/Fields/soca_fields_mod.F90

@@ -1507,13 +1507,11 @@ end function soca_genfilename
 ! ------------------------------------------------------------------------------
 !> Get the fields listed in vars, used by the interpolation.
 !!
-!! The fields that are returned 1) have halos and 2) have had the masked points
-!! removed if masked==true.
-subroutine soca_fields_to_fieldset(self, vars, afieldset, masked)
+!! The fields that are returned have halos (minus the invalid and duplicate halo points)
+subroutine soca_fields_to_fieldset(self, vars, afieldset)
   class(soca_fields),   intent(in)    :: self
   type(oops_variables), intent(in)    :: vars
   type(atlas_fieldset), intent(inout) :: afieldset
-  logical,              intent(in)    :: masked
 
   type(atlas_field) :: afield
   integer :: v, z
@@ -1528,21 +1526,6 @@ subroutine soca_fields_to_fieldset(self, vars, afieldset, masked)
     ! make sure halos are updated (remove? is redundant?)
     call field%update_halo(self%geom)
 
-    ! which mask to use
-    nullify(mask)
-    if (masked .and. field%metadata%masked) then
-      select case(field%metadata%grid)
-      case ('h')
-        mask => self%geom%mask2d
-      case ('u')
-        mask => self%geom%mask2du
-      case ('v')
-        mask => self%geom%mask2dv
-      case default
-        call abor1_ftn('incorrect grid type in soca_fields_to_fieldset()')
-      end select
-    end if
-
     ! get/create field
     if (afieldset%has_field(vars%variable(v))) then
       afield = afieldset%field(vars%variable(v))
@@ -1551,17 +1534,16 @@ subroutine soca_fields_to_fieldset(self, vars, afieldset, masked)
         name=vars%variable(v), kind=atlas_real(kind_real), levels=field%nz)
       meta = afield%metadata()
       call meta%set('interp_type', 'default')
+      if (field%metadata%masked) then
+        call meta%set('interp_source_point_mask', 'interp_mask')
+      end if
       call afieldset%add(afield)
     end if
 
     ! create and fill field
     call afield%data(real_ptr)
     do z=1,field%nz
-      if (associated(mask)) then
-        real_ptr(z,:) = pack(field%val(:,:, z), mask=mask/=0)
-      else
-        real_ptr(z,:) = pack(field%val(:,:, z), mask=.true.)
-      end if
+        real_ptr(z,:) = pack(field%val(:,:, z), mask=self%geom%valid_halo_mask)
     end do
     call afield%final()
 
@@ -1571,19 +1553,16 @@ subroutine soca_fields_to_fieldset(self, vars, afieldset, masked)
 ! ------------------------------------------------------------------------------
 !> Adjoint of get fields used by the interpolation.
 !!
-!! The fields that are input 1) have halos and 2) have had the masked points
-!! removed.
-subroutine soca_fields_to_fieldset_ad(self, vars, afieldset, masked)
+!! 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
-  logical,              intent(in) :: masked
 
   integer :: v, z
   integer :: is, ie, js, je
   type(soca_field), pointer :: field
   type(atlas_field) :: afield
-  real(kind=kind_real), pointer :: mask(:,:) => null() !< field mask
   real(kind=kind_real), pointer :: real_ptr(:,:)
   real(kind=kind_real), pointer :: tmp(:,:)
 
@@ -1597,29 +1576,10 @@ subroutine soca_fields_to_fieldset_ad(self, vars, afieldset, masked)
     call self%get(vars%variable(v), field)
     afield = afieldset%field(vars%variable(v))
 
-    ! which mask to use
-    nullify(mask)
-    if (masked .and. field%metadata%masked) then
-      select case(field%metadata%grid)
-      case ('h')
-        mask => self%geom%mask2d
-      case ('u')
-        mask => self%geom%mask2du
-      case ('v')
-        mask => self%geom%mask2dv
-      case default
-        call abor1_ftn('incorrect grid type in soca_fields_to_fieldset_ad()')
-      end select
-    end if
-
     tmp = 0.0
     call afield%data(real_ptr)
     do z=1,field%nz
-      if (associated(mask)) then
-        tmp = unpack(real_ptr(z,:), mask/=0, tmp)
-      else
-        tmp = reshape(real_ptr(z,:), shape(tmp))
-      end if
+      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
@@ -1631,25 +1591,18 @@ subroutine soca_fields_to_fieldset_ad(self, vars, afieldset, masked)
 
 ! ------------------------------------------------------------------------------
 !> Set the our values from an atlas fieldset
-subroutine soca_fields_from_fieldset(self, vars, afieldset, masked)
+subroutine soca_fields_from_fieldset(self, vars, afieldset)
   class(soca_fields), target, intent(inout) :: self
   type(oops_variables),       intent(in)    :: vars
   type(atlas_fieldset),       intent(in)    :: afieldset
-  logical,                    intent(in)    :: masked
 
-  integer :: jvar, i, jz, ngrid(2)
+  integer :: jvar, i, jz
   real(kind=kind_real), pointer :: real_ptr(:,:)
   logical :: var_found
   character(len=1024) :: fieldname
   type(soca_field), pointer :: field
   type(atlas_field) :: afield
 
-  if (masked) then
-    call abor1_ftn('soca_fields_from_fieldset() does not support masked=true')
-  end if
-
-  ngrid = (/self%geom%ied-self%geom%isd+1, self%geom%jed-self%geom%jsd+1/)
-
   ! Initialization
   call self%zeros()
 
@@ -1664,7 +1617,8 @@ subroutine soca_fields_from_fieldset(self, vars, afieldset, masked)
         ! Copy data
         call afield%data(real_ptr)
         do jz=1,field%nz
-          field%val(:,:,jz) = reshape(real_ptr(jz,:), ngrid)
+          ! 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)
         end do
 
         ! Release pointer

src/soca/Geometry/Geometry.cc

@@ -47,46 +47,6 @@ namespace soca {
       soca_geo_gridgen_f90(keyGeom_);
       return;
     }
-
-    // create kdtrees
-    int kdidx = 0;
-    for (auto grid : grids) {
-      std::vector<double> lats;
-      std::vector<double> lons;
-      size_t npoints;
-      std::vector<size_t> indx;
-
-      // kd tree with all grid points
-      this->latlon(lats, lons, true, grid, false);
-      npoints = lats.size();
-      indx.resize(npoints);
-      for (size_t jj = 0; jj < npoints; ++jj) indx[jj] = jj;
-      localTree_[kdidx++].build(lons, lats, indx);
-
-      // kd tree with only masked points
-      // NOTE: the index from the tree still refers to grid index
-      // for fields with ALL gridpoints.
-      std::vector<double> lats_masked;
-      std::vector<double> lons_masked;
-      this->latlon(lats_masked, lons_masked, true, grid, true);
-      std::vector<double> mask;
-      this->mask(mask, true, grid);
-      size_t npoints_masked = lats_masked.size();
-      if (npoints_masked > 0) {
-        // only create the tree if there are valid ocean points
-        indx.resize(npoints_masked);
-        size_t idx = 0;
-        for (size_t jj = 0; jj < npoints; ++jj) {
-          if (mask[jj] == 0.0) continue;  // land, skip
-          ASSERT(idx < npoints_masked);
-          indx[idx++] = jj;
-        }
-        ASSERT(idx == npoints_masked);
-        localTree_[kdidx++].build(lons_masked, lats_masked, indx);
-      } else {
-        kdidx++;
-      }
-    }
   }
   // -----------------------------------------------------------------------------
   Geometry::Geometry(const Geometry & other)
@@ -150,67 +110,14 @@ namespace soca {
   // -----------------------------------------------------------------------------
   void Geometry::latlon(std::vector<double> & lats, std::vector<double> & lons,
       const bool halo) const {
-    // Assume that we can get by with just using the unmasked H grid here (i.e.
-    // ignore the different staggered grids)
-    // (This method should only be called by code in OOPS used for determining the
-    //  PE that owns a specific point.... so it doesn't have to be exact)
-    latlon(lats, lons, halo, 'h', false);
-  }
-  // -----------------------------------------------------------------------------
-  void Geometry::latlon(
-        std::vector<double> & lats, std::vector<double> & lons, const bool halo,
-        const char grid, const bool masked) const {
     // get the number of gridpoints
     int gridSize;
-    soca_geo_gridsize_f90(keyGeom_, grid, masked, halo, gridSize);
+    soca_geo_gridsize_f90(keyGeom_, halo, gridSize);
 
     // get the lat/lon of those gridpoints
     lats.resize(gridSize);
     lons.resize(gridSize);
-    soca_geo_gridlatlon_f90(keyGeom_, grid, masked, halo, gridSize,
-      lats.data(), lons.data());
-  }
-  // -----------------------------------------------------------------------------
-  void Geometry::mask(
-      std::vector<double> & mask, const bool halo, const char grid) const {
-    // get the number of gridpoints
-    int gridSize;
-    soca_geo_gridsize_f90(keyGeom_, grid, false, halo, gridSize);
-    mask.resize(gridSize);
-
-    // get the mask value of those gridpoints
-    soca_geo_gridmask_f90(keyGeom_, grid, halo, gridSize, mask.data());
-  }
-
-  // -----------------------------------------------------------------------------
-  atlas::util::KDTree<size_t>::ValueList Geometry::closestPoints(
-    const double lat, const double lon, const int npoints,
-    const char grid, const bool masked) const {
-
-    // determine which kdtree to use
-    // TODO(travis) make sure not -1
-    int kdidx = -1;
-    for (int j=0; j < grids.size(); j++) {
-      if (grids[j] == grid) kdidx = j*2;
-    }
-    ASSERT(kdidx >= 0);
-    if (masked) kdidx += 1;
-    ASSERT(kdidx < 6);
-
-    atlas::PointLonLat obsloc(lon, lat);
-    obsloc.normalise();
-    atlas::util::KDTree<size_t>::ValueList neighbours =
-      localTree_[kdidx].closestPoints(obsloc, npoints);
-
-    return neighbours;
-    }
-  // -----------------------------------------------------------------------------
-  // Get the properties of the grid for a given variable (masked and u/v/h grid)
-  void Geometry::getVarGrid(const std::string &var, char & grid, bool & masked) const {
-      oops::Variables vars;
-      vars.push_back(var);
-
-      soca_geo_getvargrid_f90(keyGeom_, vars, grid, masked);
+    soca_geo_gridlatlon_f90(keyGeom_, halo, gridSize, lats.data(), lons.data());
   }
 
 }  // namespace soca

src/soca/Geometry/Geometry.h

@@ -75,13 +75,6 @@ namespace soca {
       atlas::FieldSet & extraFields() {return extraFields_;}
 
       void latlon(std::vector<double> &, std::vector<double> &, const bool) const;
-      void latlon(std::vector<double> &, std::vector<double> &, const bool,
-                  const char, const bool) const;
-      atlas::util::KDTree<size_t>::ValueList closestPoints(
-        const double, const double, const int,
-        const char, const bool) const;
-
-      void getVarGrid(const std::string &, char &, bool &) const;
 
    private:
       void mask(std::vector<double> &, const bool, const char) const;
@@ -94,8 +87,6 @@ namespace soca {
       atlas::FunctionSpace functionSpace_;
       atlas::FunctionSpace functionSpaceIncHalo_;
       atlas::FieldSet extraFields_;
-
-      atlas::util::IndexKDTree localTree_[6];
   };
   // -----------------------------------------------------------------------------
 

src/soca/Geometry/GeometryFortran.h

@@ -49,14 +49,8 @@ namespace soca {
                                     const size_t &, size_t[]);
     void soca_geo_iterator_dimension_f90(const F90geom &, int &);
 
-    void soca_geo_gridsize_f90(const F90geom &, const char &, const bool &,
-      const bool &, int &);
-    void soca_geo_gridlatlon_f90(const F90geom &, const char &, const bool &,
-      const bool &, const int &, double[], double[]);
-    void soca_geo_gridmask_f90(const F90geom&, const char &, const bool &,
-      const int &, double[]);
-    void soca_geo_getvargrid_f90(
-      const F90geom &, const oops::Variables &, char &, bool &);
+    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_