Add capability to extract slices from Fields in Fortran

src/atlas/field/FieldSet.cc

@@ -11,6 +11,7 @@
 #include <sstream>
 
 #include "atlas/field/Field.h"
+#include "atlas/field/detail/FieldInterface.h"
 #include "atlas/field/FieldSet.h"
 #include "atlas/grid/Grid.h"
 #include "atlas/runtime/Exception.h"
@@ -86,6 +87,39 @@ std::vector<std::string> FieldSetImpl::field_names() const {
 // C wrapper interfaces to C++ routines
 extern "C" {
 
+void atlas__FieldSet__data_int_specf(FieldSetImpl* This, char* name, int*& data, int& rank, int*& shapef, int*& stridesf) {
+    atlas__Field__data_int_specf(This->field(name).get(), data, rank, shapef, stridesf);
+}
+
+void atlas__FieldSet__data_long_specf(FieldSetImpl* This, char* name, long*& data, int& rank, int*& shapef, int*& stridesf) {
+    atlas__Field__data_long_specf(This->field(name).get(), data, rank, shapef, stridesf);
+}
+
+void atlas__FieldSet__data_float_specf(FieldSetImpl* This, char* name, float*& data, int& rank, int*& shapef, int*& stridesf) {
+    atlas__Field__data_float_specf(This->field(name).get(), data, rank, shapef, stridesf);
+}
+
+void atlas__FieldSet__data_double_specf(FieldSetImpl* This, char* name, double*& data, int& rank, int*& shapef, int*& stridesf) {
+    atlas__Field__data_double_specf(This->field(name).get(), data, rank, shapef, stridesf);
+}
+
+void atlas__FieldSet__data_int_specf_by_idx(FieldSetImpl* This, int& idx, int*& data, int& rank, int*& shapef, int*& stridesf) {
+    atlas__Field__data_int_specf(This->operator[](idx).get(), data, rank, shapef, stridesf);
+}
+
+void atlas__FieldSet__data_long_specf_by_idx(FieldSetImpl* This, int& idx, long*& data, int& rank, int*& shapef, int*& stridesf) {
+    atlas__Field__data_long_specf(This->operator[](idx).get(), data, rank, shapef, stridesf);
+}
+
+void atlas__FieldSet__data_float_specf_by_idx(FieldSetImpl* This, int& idx, float*& data, int& rank, int*& shapef, int*& stridesf) {
+    atlas__Field__data_float_specf(This->operator[](idx).get(), data, rank, shapef, stridesf);
+}
+
+void atlas__FieldSet__data_double_specf_by_idx(FieldSetImpl* This, int& idx, double*& data, int& rank, int*& shapef, int*& stridesf) {
+    atlas__Field__data_double_specf(This->operator[](idx).get(), data, rank, shapef, stridesf);
+}
+
+
 FieldSetImpl* atlas__FieldSet__new(char* name) {
     FieldSetImpl* fset = new FieldSetImpl(std::string(name));
     fset->name()       = name;

src/atlas/field/FieldSet.h

@@ -24,6 +24,7 @@
 
 #include "eckit/deprecated.h"
 
+#include "atlas/array_fwd.h"
 #include "atlas/field/Field.h"
 #include "atlas/library/config.h"
 #include "atlas/runtime/Exception.h"
@@ -139,6 +140,22 @@ const char* atlas__FieldSet__name(FieldSetImpl* This);
 idx_t atlas__FieldSet__size(const FieldSetImpl* This);
 FieldImpl* atlas__FieldSet__field_by_name(FieldSetImpl* This, char* name);
 FieldImpl* atlas__FieldSet__field_by_idx(FieldSetImpl* This, idx_t idx);
+void atlas__FieldSet__data_int_specf(FieldSetImpl* This, char* name, int*& field_data, int& rank, int*& field_shapef,
+                                  int*& field_stridesf);
+void atlas__FieldSet__data_long_specf(FieldSetImpl* This, char* name, long*& field_data, int& rank, int*& field_shapef,
+                                   int*& field_stridesf);
+void atlas__FieldSet__data_float_specf(FieldSetImpl* This, char* name, float*& field_data, int& rank, int*& field_shapef,
+                                    int*& field_stridesf);
+void atlas__FieldSet__data_double_specf(FieldSetImpl* This, char* name, double*& field_data, int& rank, int*& field_shapef,
+                                     int*& field_stridesf);
+void atlas__FieldSet__data_int_specf_by_idx(FieldSetImpl* This, int& idx, int*& field_data, int& rank, int*& field_shapef,
+                                  int*& field_stridesf);
+void atlas__FieldSet__data_long_specf_by_idx(FieldSetImpl* This, int& idx, long*& field_data, int& rank, int*& field_shapef,
+                                   int*& field_stridesf);
+void atlas__FieldSet__data_float_specf_by_idx(FieldSetImpl* This, int& idx, float*& field_data, int& rank, int*& field_shapef,
+                                    int*& field_stridesf);
+void atlas__FieldSet__data_double_specf_by_idx(FieldSetImpl* This, int& idx, double*& field_data, int& rank, int*& field_shapef,
+                                     int*& field_stridesf);
 void atlas__FieldSet__set_dirty(FieldSetImpl* This, int value);
 void atlas__FieldSet__halo_exchange(FieldSetImpl* This, int on_device);
 }

src/atlas/functionspace/BlockStructuredColumns.h

@@ -42,8 +42,6 @@ class BlockStructuredColumns : public FunctionSpace {
     bool valid() const { return functionspace_; }
 
     idx_t size() const { return functionspace_->size(); }
-//    idx_t sizeOwned() const { return functionspace_->sizeOwned(); }
-//    idx_t sizeHalo() const { return functionspace_->sizeHalo(); }
     idx_t levels() const { return functionspace_->levels(); }
 
     const Vertical& vertical() const { return functionspace_->vertical(); }
@@ -54,10 +52,6 @@ class BlockStructuredColumns : public FunctionSpace {
     std::string checksum(const Field&) const;
 
     idx_t index(idx_t blk, idx_t rof) const { return functionspace_->index(blk, rof); }
-//    idx_t i_begin(idx_t j) const { return functionspace_->i_begin(j); }
-//    idx_t i_end(idx_t j) const { return functionspace_->i_end(j); }
-//    idx_t j_begin() const { return functionspace_->j_begin(); }
-//    idx_t j_end() const { return functionspace_->j_end(); }
     idx_t k_begin() const { return functionspace_->k_begin(); }
     idx_t k_end() const { return functionspace_->k_end(); }
     idx_t nproma() const { return functionspace_->nproma(); }

src/atlas/functionspace/detail/BlockStructuredColumns.cc

@@ -156,7 +156,6 @@ void rev_block_copy(const Field loc, Field sloc, const functionspace::detail::Bl
     }
 }
 
-
 void transpose_nonblocked_to_blocked(const Field& nonblocked, Field& blocked, const functionspace::detail::BlockStructuredColumns& fs) {
     auto kind = nonblocked.datatype().kind();
     if (kind == array::DataType::kind<int>()) {

src/atlas/functionspace/detail/StructuredColumns.h

@@ -317,6 +317,7 @@ class StructuredColumns : public FunctionSpaceImpl {
     friend struct BlockStructuredColumnsFortranAccess;
     Map2to1 ij2gp_;
 
+    friend class BlockStructuredColumns;
     void setup(const grid::Distribution& distribution, const eckit::Configuration& config);
 
     friend class BlockStructuredColumns;

src/atlas_f/CMakeLists.txt

@@ -175,7 +175,7 @@ ecbuild_add_library( TARGET atlas_f
         functionspace/atlas_functionspace_BlockStructuredColumns_module.F90
         functionspace/atlas_functionspace_Spectral_module.F90
         functionspace/atlas_functionspace_PointCloud_module.F90
-        field/atlas_FieldSet_module.F90
+        field/atlas_FieldSet_module.fypp
         field/atlas_State_module.F90
         field/atlas_Field_module.fypp
         grid/atlas_Grid_module.F90

src/atlas_f/field/atlas_FieldSet_module.F90 → src/atlas_f/field/atlas_FieldSet_module.fypp

@@ -8,9 +8,13 @@
 
 #include "atlas/atlas_f.h"
 
+#:include "atlas/atlas_f.fypp"
+#:include "internals/atlas_generics.fypp"
+
 module atlas_FieldSet_module
 
 use fckit_owned_object_module, only: fckit_owned_object
+use atlas_field_module, only: atlas_field, array_c_to_f
 use atlas_kinds_module, only : ATLAS_KIND_IDX
 
 implicit none
@@ -53,10 +57,33 @@ module atlas_FieldSet_module
   procedure, public :: set_dirty
   procedure, public :: halo_exchange
 
+#:for rank in ranks
+#:for dtype in dtypes
+  procedure, private :: access_data_${dtype}$_r${rank}$_by_name
+  procedure, private :: access_data_${dtype}$_r${rank}$_by_idx
+  procedure, private :: access_data_${dtype}$_r${rank}$_slice_by_name
+  procedure, private :: access_data_${dtype}$_r${rank}$_slice_by_idx
+#:endfor
+#:endfor
+
+  generic, public :: data => &
+#:for rank in ranks
+#:for dtype in dtypes
+      & access_data_${dtype}$_r${rank}$_by_name, &
+      & access_data_${dtype}$_r${rank}$_by_idx, &
+      & access_data_${dtype}$_r${rank}$_slice_by_name, &
+      & access_data_${dtype}$_r${rank}$_slice_by_idx, &
+#:endfor
+#:endfor
+      & dummy
+
+  procedure, private :: dummy
+
 #if FCKIT_FINAL_NOT_INHERITING
   final :: atlas_FieldSet__final_auto
 #endif
   procedure, public  :: has_field => has ! deprecated !
+
 END TYPE atlas_FieldSet
 !------------------------------------------------------------------------------
 
@@ -74,6 +101,86 @@ module atlas_FieldSet_module
 ! -----------------------------------------------------------------------------
 ! FieldSet routines
 
+#:for rank  in ranks
+#:for dtype,ftype,ctype in types
+subroutine access_data_${dtype}$_r${rank}$_by_name(this, name, field)
+  use fckit_c_interop_module, only: c_str
+  use atlas_fieldset_c_binding
+  use, intrinsic :: iso_c_binding, only : c_ptr, c_int, c_long, c_float, c_double
+  class(atlas_FieldSet), intent(in) :: this
+  character(len=*), intent(in) :: name
+  ${ftype}$, pointer, intent(inout) :: field(${dim[rank]}$)
+  type(c_ptr) :: field_cptr
+  type(c_ptr) :: shape_cptr
+  type(c_ptr) :: strides_cptr
+  integer(c_int) :: rank
+  call atlas__FieldSet__data_${ctype}$_specf(this%CPTR_PGIBUG_A, c_str(name), field_cptr, rank, shape_cptr, strides_cptr)
+  call array_c_to_f(field_cptr, rank, shape_cptr, strides_cptr, field)
+end subroutine
+subroutine access_data_${dtype}$_r${rank}$_by_idx(this, idx, field)
+  use fckit_c_interop_module, only: c_str
+  use atlas_fieldset_c_binding
+  use, intrinsic :: iso_c_binding, only : c_ptr, c_int, c_long, c_float, c_double
+  class(atlas_FieldSet), intent(in) :: this
+  integer, intent(in) :: idx
+  ${ftype}$, pointer, intent(inout) :: field(${dim[rank]}$)
+  type(c_ptr) :: field_cptr
+  type(c_ptr) :: shape_cptr
+  type(c_ptr) :: strides_cptr
+  integer(c_int) :: rank
+  call atlas__FieldSet__data_${ctype}$_specf_by_idx(this%CPTR_PGIBUG_A, idx-1, field_cptr, rank, shape_cptr, strides_cptr)
+  call array_c_to_f(field_cptr, rank, shape_cptr, strides_cptr, field)
+end subroutine
+subroutine access_data_${dtype}$_r${rank}$_slice_by_name(this, name, slice, iblk)
+  use fckit_c_interop_module, only: c_str
+  use atlas_fieldset_c_binding
+  use, intrinsic :: iso_c_binding, only : c_ptr, c_int, c_long, c_float, c_double
+  class(atlas_FieldSet), intent(in) :: this
+  character(len=*), intent(in) :: name
+#:if rank > 1
+  ${ftype}$, pointer, intent(inout) :: slice(${dimr[rank]}$)
+#:else
+  ${ftype}$, pointer, intent(inout) :: slice
+#:endif
+  integer, intent(in) :: iblk
+  ${ftype}$, pointer :: field(${dim[rank]}$)
+  call access_data_${dtype}$_r${rank}$_by_name(this, c_str(name), field)
+#:if rank > 1
+  slice => field(${dimr[rank]}$, iblk)
+#:else
+  slice => field(iblk)
+#:endif 
+end subroutine
+subroutine access_data_${dtype}$_r${rank}$_slice_by_idx(this, idx, slice, iblk)
+  use fckit_c_interop_module, only: c_str
+  use atlas_fieldset_c_binding
+  use, intrinsic :: iso_c_binding, only : c_ptr, c_int, c_long, c_float, c_double
+  class(atlas_FieldSet), intent(in) :: this
+  integer, intent(in) :: idx
+#:if rank > 1
+  ${ftype}$, pointer, intent(inout) :: slice(${dimr[rank]}$)
+#:else
+  ${ftype}$, pointer, intent(inout) :: slice
+#:endif
+  integer, intent(in) :: iblk
+  ${ftype}$, pointer :: field(${dim[rank]}$)
+  call access_data_${dtype}$_r${rank}$_by_idx(this, idx, field)
+#:if rank > 1
+  slice => field(${dimr[rank]}$, iblk)
+#:else
+  slice => field(iblk)
+#:endif 
+end subroutine
+!-------------------------------------------------------------------------------
+
+#:endfor
+#:endfor
+subroutine dummy(this)
+  use atlas_fieldset_c_binding
+  class(atlas_FieldSet), intent(in) :: this
+  FCKIT_SUPPRESS_UNUSED(this)
+end subroutine
+
 function atlas_FieldSet__cptr(cptr) result(fieldset)
   use, intrinsic :: iso_c_binding, only: c_ptr
   type(atlas_FieldSet) :: fieldset

src/atlas_f/field/atlas_Field_module.fypp

@@ -23,6 +23,7 @@ public :: atlas_real
 public :: atlas_integer
 public :: atlas_logical
 public :: atlas_data_type
+public :: array_c_to_f
 
 private
 
@@ -75,6 +76,7 @@ contains
 #:for dtype in dtypes
   procedure, private :: access_data_${dtype}$_r${rank}$
   procedure, private :: access_data_${dtype}$_r${rank}$_shape
+  procedure, private :: access_data_${dtype}$_r${rank}$_slice
 #:endfor
 #:endfor
 
@@ -83,6 +85,7 @@ contains
 #:for dtype in dtypes
       & access_data_${dtype}$_r${rank}$, &
       & access_data_${dtype}$_r${rank}$_shape, &
+      & access_data_${dtype}$_r${rank}$_slice, &
 #:endfor
 #:endfor
       & dummy
@@ -135,7 +138,6 @@ interface array_c_to_f
 end interface
 !-------------------------------------------------------------------------------
 
-
 private :: fckit_owned_object
 private :: atlas_Config
 
@@ -161,6 +163,7 @@ subroutine array_c_to_f_${dtype}$_r${rank}$(array_cptr,rank,shape_cptr,strides_c
   integer :: j
 
   if( rank /= ${rank}$ ) then
+    write(0,*) rank, "!=", ${rank}$
     @:ATLAS_ABORT("Rank mismatch")
   endif
 
@@ -199,6 +202,29 @@ subroutine access_data_${dtype}$_r${rank}$(this, field)
   call atlas__Field__data_${ctype}$_specf(this%CPTR_PGIBUG_A, field_cptr, rank, shape_cptr, strides_cptr)
   call array_c_to_f(field_cptr,rank,shape_cptr,strides_cptr, field)
 end subroutine
+subroutine access_data_${dtype}$_r${rank}$_slice(this, slice, iblk)
+  use atlas_field_c_binding
+  use, intrinsic :: iso_c_binding, only : c_ptr, c_int, c_long, c_float, c_double
+  class(atlas_Field), intent(in) :: this
+  ${ftype}$, pointer :: field(${dim[rank]}$)
+#:if rank > 1
+  ${ftype}$, pointer, intent(inout) :: slice(${dimr[rank]}$)
+#:else
+  ${ftype}$, pointer, intent(inout) :: slice
+#:endif
+  integer, intent(in) :: iblk
+  type(c_ptr) :: field_cptr
+  type(c_ptr) :: shape_cptr
+  type(c_ptr) :: strides_cptr
+  integer(c_int) :: rank
+  call atlas__Field__data_${ctype}$_specf(this%CPTR_PGIBUG_A, field_cptr, rank, shape_cptr, strides_cptr)
+  call array_c_to_f(field_cptr, rank, shape_cptr, strides_cptr, field)
+#:if rank > 1
+  slice => field(${dimr[rank]}$, iblk)
+#:else
+  slice => field(iblk)
+#:endif
+end subroutine
 
 !-------------------------------------------------------------------------------
 

src/atlas_f/internals/atlas_generics.fypp

@@ -10,6 +10,7 @@
 
 #:set ranks  = [1,2,3,4]
 #:set dim    = ['',':',':,:',':,:,:',':,:,:,:',':,:,:,:,:']
+#:set dimr   = ['','',':',':,:',':,:,:',':,:,:,:']
 
 #:set ftypes = ['integer(c_int)','integer(c_long)','real(c_float)','real(c_double)', 'logical']
 #:set ctypes = ['int','long','float','double', 'int']