diff --git a/doc/CMakeLists.txt b/doc/CMakeLists.txt
index 6c7fa7d3a..05c754c4b 100644
--- a/doc/CMakeLists.txt
+++ b/doc/CMakeLists.txt
@@ -6,6 +6,10 @@
 # granted to it by virtue of its status as an intergovernmental organisation nor
 # does it submit to any jurisdiction.
 
+if(atlas_HAVE_FORTRAN)
+  add_subdirectory(example-fortran)
+endif()
+
 if(atlas_HAVE_DOCS)
 
 add_subdirectory(user-guide)
diff --git a/doc/example-fortran/CMakeLists.txt b/doc/example-fortran/CMakeLists.txt
new file mode 100644
index 000000000..564a289f0
--- /dev/null
+++ b/doc/example-fortran/CMakeLists.txt
@@ -0,0 +1,10 @@
+# (C) Copyright 2024- ECMWF.
+#
+# 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.
+# In applying this licence, ECMWF does not waive the privileges and immunities
+# granted to it by virtue of its status as an intergovernmental organisation nor
+# does it submit to any jurisdiction.
+
+add_subdirectory(structured-grid-stencils)
+
diff --git a/doc/example-fortran/structured-grid-stencils/CMakeLists.txt b/doc/example-fortran/structured-grid-stencils/CMakeLists.txt
new file mode 100644
index 000000000..c3a221662
--- /dev/null
+++ b/doc/example-fortran/structured-grid-stencils/CMakeLists.txt
@@ -0,0 +1,12 @@
+# (C) Copyright 2024- ECMWF.
+#
+# 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.
+# In applying this licence, ECMWF does not waive the privileges and immunities
+# granted to it by virtue of its status as an intergovernmental organisation nor
+# does it submit to any jurisdiction.
+
+ecbuild_add_executable( TARGET example-fortran-structured-grid-stencils
+                        OUTPUT_NAME program NOINSTALL
+                        SOURCES program.F90
+                        LIBS atlas_f )
diff --git a/doc/example-fortran/structured-grid-stencils/program.F90 b/doc/example-fortran/structured-grid-stencils/program.F90
new file mode 100644
index 000000000..179fcd51a
--- /dev/null
+++ b/doc/example-fortran/structured-grid-stencils/program.F90
@@ -0,0 +1,208 @@
+! (C) Copyright 2023 ECMWF.
+!
+! 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.
+! In applying this licence, ECMWF does not waive the privileges and immunities
+! granted to it by virtue of its status as an intergovernmental organisation nor
+! does it submit to any jurisdiction.
+
+! ------------------------------------------------------------------------------------------------
+! Example program on how to access and use StructuredColumns functionspace
+! Also includes computation of stencils around a given coordinate
+! ------------------------------------------------------------------------------------------------
+
+subroutine run
+
+use, intrinsic :: iso_fortran_env, only : real64, real32
+use atlas_module, only :   &
+  & atlas_StructuredGrid  ,&
+  & atlas_functionspace_StructuredColumns ,&
+  & atlas_Field, &
+  & ATLAS_KIND_IDX, &
+  & atlas_StructuredGrid_ComputeNorth ,&
+  & atlas_StructuredGrid_ComputeWest ,&
+  & atlas_StructuredGrid_ComputeStencil ,&
+  & atlas_StructuredGrid_Stencil, &
+  & atlas_real
+
+use fckit_mpi_module, only : fckit_mpi
+
+! ------------------------------------------------
+implicit none
+
+character(len=:), allocatable :: gridname
+type(atlas_StructuredGrid) :: grid
+type(atlas_functionspace_StructuredColumns) :: functionspace
+real(real64) :: zlon, zlat
+integer(ATLAS_KIND_IDX) :: jlon,jlat, jgp
+real(real64), pointer :: xy(:,:)
+integer :: log_rank
+
+! ------------------------------------------------
+
+log_rank = 1 ! First partition has rank=0 !!!
+if (fckit_mpi%size() == 1) then
+  log_rank = 0
+endif
+
+gridname = "O320"
+
+! Create grid and write some properties
+grid          = atlas_StructuredGrid(gridname)
+
+if (fckit_mpi%rank() == log_rank) then
+    write(0,*) "grid%name()  : ", grid%name()
+    write(0,*) "grid%size()  : ", grid%size()
+    write(0,*) "grid%ny()    : ", grid%ny()
+    write(0,*) "grid%nxmax() : ", grid%nxmax()
+endif
+
+! Create functionspace of grid, distributed across MPI tasks with default partitioner, and with halo of 2 wide
+!   The default partitioner is "equal_regions"
+functionspace = atlas_functionspace_StructuredColumns(grid, halo=2)
+
+! Latitude bounds without and with halos of this partition
+if (fckit_mpi%rank() == log_rank) then
+    write(0,*) "functionspace%size_owned() : ", functionspace%size_owned()
+    write(0,*) "functionspace%size()       : ", functionspace%size()
+    write(0,'(A,I0,A,I0)') " halo points between index ", functionspace%size_owned()+1 , " and " , functionspace%size()
+    write(0,*) "functionspace%j_begin(),      functionspace%j_end()      : ", &
+      & functionspace%j_begin(), functionspace%j_end()
+    write(0,*) "functionspace%j_begin_halo(), functionspace%j_end_halo() : ", &
+      & functionspace%j_begin_halo(), functionspace%j_end_halo()
+endif
+
+! Longitude bounds without and with halos of this partition for the first latitude
+jlat = functionspace%j_begin()
+if (fckit_mpi%rank() == log_rank) then
+  write(0,*) "functionspace%i_begin(jlat),      functionspace%i_end(jlat)      : ", functionspace%i_begin(jlat),     &
+   & functionspace%i_end(jlat)
+  write(0,*) "functionspace%i_begin_halo(jlat), functionspace%i_end_halo(jlat) : ", functionspace%i_begin_halo(jlat),&
+   & functionspace%i_end_halo(jlat)
+endif
+
+! Access i,j indices of grid points
+block
+  type(atlas_Field) :: field_index_i, field_index_j
+  integer(ATLAS_KIND_IDX), pointer :: index_i(:), index_j(:)
+  field_index_i = functionspace%index_i()
+  field_index_j = functionspace%index_j()
+  call field_index_i%data(index_i)
+  call field_index_j%data(index_j)
+  if (fckit_mpi%rank() == log_rank) then
+    write(0,*) "i,j of first partition point :", index_i(1), index_j(1)
+  endif
+  call field_index_i%final()
+  call field_index_j%final()
+end block
+
+! Access to xy coordinates
+block
+  type(atlas_Field) :: field_xy
+  field_xy = functionspace%xy()
+  call field_xy%data(xy)
+  call field_xy%final()
+end block
+
+! Creating a horizontal field and perform halo exchange
+block
+  type(atlas_Field) :: field
+  real(real32), pointer :: view(:)
+  field = functionspace%create_field(name="myfield", kind=atlas_real(real32))
+  call field%data(view)
+  if (fckit_mpi%rank() == log_rank) then
+    write(0,*) "shape( horizontal field )", shape(view)
+  endif
+  view(1:functionspace%size_owned()) = 1.
+  view(functionspace%size_owned()+1:functionspace%size()) = 0
+  call field%set_dirty() ! This marks that halos are in "dirty" state
+  call field%halo_exchange() ! Only exchanges halos when halos are marked as "dirty"
+  if (fckit_mpi%rank() == log_rank) then
+    write(0,*) "halo exhange success : ", all(view == 1.)
+  endif
+  call field%final()
+end block
+
+! Creating a horizontal/vertical field
+block
+  type(atlas_Field) :: field
+  real(real32), pointer :: view(:,:)
+  field = functionspace%create_field(name="myfield", kind=atlas_real(real32), levels=10)
+  call field%data(view)
+  if (fckit_mpi%rank() == log_rank) then
+    write(0,*) "shape( horizontal/vertical field )", shape(view)
+  endif
+  call field%final()
+end block
+
+! Set a coordinate somewhere south-east of the first grid point, but still within this partition
+jgp = 1
+zlon = xy(1,jgp) + 0.1
+zlat = xy(2,jgp) - 0.1
+
+! Compute nearest points to the north-west of a coordinate
+block
+  type(atlas_StructuredGrid_ComputeNorth) :: compute_north
+  type(atlas_StructuredGrid_ComputeWest) :: compute_west
+  compute_north = atlas_StructuredGrid_ComputeNorth(grid, halo=2)
+  compute_west  = atlas_StructuredGrid_ComputeWest(grid,  halo=2)
+  
+  jlat = compute_north%execute(zlat)
+  jlon = compute_west%execute(zlon, jlat)
+  jgp = functionspace%index(jlon,jlat) ! gridpoint index in 1-dimensional array
+  
+  if (fckit_mpi%rank() == log_rank) then
+  write(0,'(A,F6.2,A,I0)') "compute_north%execute(y=",zlat,") : ", jlat
+  write(0,'(A,F6.2,A,I0,A,I0)') "compute_west%execute(x=",zlon,", j=", jlat,") : ", jlon
+  write(0,'(A,I0,A,F6.2,A,F6.2,A)') "gridpoint north-west: jpg=",jgp,"    (lon,lat)=(", xy(1,jgp), ",", xy(2,jgp), ")"
+  endif
+  call compute_west%final()
+  call compute_north%final()
+end block
+
+! Stencil computer of 4x4 stencil around a coordinate
+block
+  type(atlas_StructuredGrid_ComputeStencil) :: compute_stencil
+  type(atlas_StructuredGrid_Stencil) :: stencil
+  call compute_stencil%setup(grid, stencil_width=4)
+  call compute_stencil%execute(zlon,zlat,stencil)
+  if (fckit_mpi%rank() == log_rank) then
+    write(0,'(A,A,dt)') 'stencil:', new_line('a'), stencil
+  endif
+  if (fckit_mpi%rank() == log_rank) then
+    write(0,'(A,A,dt)') 'stencil gridpoints:'
+    do jlat=1,stencil%width
+      do jlon=1,stencil%width
+        write(0,'(I8)',advance='no') functionspace%index(stencil%i(jlon,jlat),stencil%j(jlat))
+      enddo
+      write(0,'(A)') ""
+    enddo
+
+    write(0,'(A,A,dt)') 'stencil coordinates:'
+    do jlat=1,stencil%width
+      jgp = functionspace%index(stencil%i(1,jlat),stencil%j(jlat))
+      write(0,'(A, F6.2,A)', advance='no') " lat : ", xy(2,jgp), "  ---  lon : "
+      do jlon=1,stencil%width
+        jgp = functionspace%index(stencil%i(jlon,jlat),stencil%j(jlat))
+        write(0,'(F8.2)',advance='no') xy(1,jgp)
+      enddo
+      write(0,'(A)') ""
+    enddo
+
+  endif
+  call compute_stencil%final()
+end block
+
+call functionspace%final()
+call grid%final()
+end subroutine
+
+! ------------------------------------------------------------------------------------------------
+
+program main
+use atlas_module, only : atlas_library
+implicit none
+call atlas_library%initialize()
+call run()
+call atlas_library%finalize()
+end program
\ No newline at end of file
diff --git a/src/atlas/CMakeLists.txt b/src/atlas/CMakeLists.txt
index 3636d74d6..1f8306cf9 100644
--- a/src/atlas/CMakeLists.txt
+++ b/src/atlas/CMakeLists.txt
@@ -173,6 +173,8 @@ grid/detail/grid/RegionalVariableResolution.h
 grid/detail/grid/RegionalVariableResolution.cc
 grid/detail/grid/Healpix.h
 grid/detail/grid/Healpix.cc
+grid/detail/grid/StencilComputerInterface.h
+grid/detail/grid/StencilComputerInterface.cc
 
 grid/detail/tiles/Tiles.cc
 grid/detail/tiles/Tiles.h
diff --git a/src/atlas/grid/detail/grid/StencilComputerInterface.cc b/src/atlas/grid/detail/grid/StencilComputerInterface.cc
new file mode 100644
index 000000000..4e1df4765
--- /dev/null
+++ b/src/atlas/grid/detail/grid/StencilComputerInterface.cc
@@ -0,0 +1,48 @@
+/*
+ * (C) Copyright 2023 ECMWF.
+ *
+ * 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.
+ * In applying this licence, ECMWF does not waive the privileges and immunities
+ * granted to it by virtue of its status as an intergovernmental organisation
+ * nor does it submit to any jurisdiction.
+ */
+
+#include "StencilComputerInterface.h"
+
+extern "C" {
+    atlas::grid::ComputeNorth* atlas__grid__ComputeNorth__new(const atlas::StructuredGrid::Implementation* grid, int halo) {
+        return new atlas::grid::ComputeNorth{atlas::StructuredGrid{grid}, halo};
+    }
+
+    void atlas__grid__ComputeNorth__delete(atlas::grid::ComputeNorth* This) {
+        delete This;
+    }
+
+    atlas::idx_t atlas__grid__ComputeNorth__execute_real32(const atlas::grid::ComputeNorth* This, float y) {
+        return This->operator()(y);
+    }
+
+    atlas::idx_t atlas__grid__ComputeNorth__execute_real64(const atlas::grid::ComputeNorth* This, double y) {
+        return This->operator()(y);
+    }
+
+
+    atlas::grid::ComputeWest* atlas__grid__ComputeWest__new(const atlas::StructuredGrid::Implementation* grid, int halo) {
+        return new atlas::grid::ComputeWest{atlas::StructuredGrid{grid}, halo};
+    }
+
+    void atlas__grid__ComputeWest__delete(atlas::grid::ComputeWest* This) {
+        delete This;
+    }
+
+    atlas::idx_t atlas__grid__ComputeWest__execute_real32(const atlas::grid::ComputeWest* This, float x, atlas::idx_t j) {
+        return This->operator()(x, j);
+    }
+
+    atlas::idx_t atlas__grid__ComputeWest__execute_real64(const atlas::grid::ComputeWest* This, double x, atlas::idx_t j) {
+        return This->operator()(x, j);
+    }
+
+
+}
diff --git a/src/atlas/grid/detail/grid/StencilComputerInterface.h b/src/atlas/grid/detail/grid/StencilComputerInterface.h
new file mode 100644
index 000000000..87e4172d2
--- /dev/null
+++ b/src/atlas/grid/detail/grid/StencilComputerInterface.h
@@ -0,0 +1,27 @@
+/*
+ * (C) Copyright 2023 ECMWF.
+ *
+ * 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.
+ * In applying this licence, ECMWF does not waive the privileges and immunities
+ * granted to it by virtue of its status as an intergovernmental organisation
+ * nor does it submit to any jurisdiction.
+ */
+
+#pragma once
+
+#include "atlas/grid/StencilComputer.h"
+#include "atlas/grid/StructuredGrid.h"
+
+extern "C" {
+    atlas::grid::ComputeNorth* atlas__grid__ComputeNorth__new(const atlas::StructuredGrid::Implementation* grid, int halo);
+    void atlas__grid__ComputeNorth__delete(atlas::grid::ComputeNorth* This);
+    atlas::idx_t atlas__grid__ComputeNorth__execute_real32(const atlas::grid::ComputeNorth* This, float y);
+    atlas::idx_t atlas__grid__ComputeNorth__execute_real64(const atlas::grid::ComputeNorth* This, double y);
+
+    atlas::grid::ComputeWest* atlas__grid__ComputeWest__new(const atlas::StructuredGrid::Implementation* grid, int halo);
+    void atlas__grid__ComputeWest__delete(atlas::grid::ComputeWest* This);
+    atlas::idx_t atlas__grid__ComputeWest__execute_real32(const atlas::grid::ComputeWest* This, float x, atlas::idx_t j);
+    atlas::idx_t atlas__grid__ComputeWest__execute_real64(const atlas::grid::ComputeWest* This, double x, atlas::idx_t j);
+
+}
diff --git a/src/atlas_f/CMakeLists.txt b/src/atlas_f/CMakeLists.txt
index 3c3d396e0..bcfd373e8 100644
--- a/src/atlas_f/CMakeLists.txt
+++ b/src/atlas_f/CMakeLists.txt
@@ -49,7 +49,7 @@ function(generate_fortran_bindings output filename)
     OUTPUT ${outfile}
     COMMAND ${PYTHON_EXECUTABLE} ${PROJECT_SOURCE_DIR}/tools/c2f.py ${CMAKE_CURRENT_SOURCE_DIR}/${filename}
             -o ${outfile} -m ${_PAR_MODULE}
-            -t '{"idx_t":"${F_IDX}","gidx_t":"${F_GIDX}"}'
+            -t '{"idx_t":"${F_IDX}","gidx_t":"${F_GIDX}","atlas::idx_t":"${F_IDX}","atlas::gidx_t":"${F_GIDX}"}'
     DEPENDS ${filename} )
   set_source_files_properties(${outfile} PROPERTIES GENERATED TRUE)
 endfunction()
@@ -137,6 +137,10 @@ generate_fortran_bindings(FORTRAN_BINDINGS ../atlas/redistribution/detail/Redist
     MODULE    atlas_redistribution_c_binding
     OUTPUT    redistribution_c_binding.f90)
 
+generate_fortran_bindings(FORTRAN_BINDINGS ../atlas/grid/detail/grid/StencilComputerInterface.h
+    MODULE    atlas_grid_StencilComputer_c_binding
+    OUTPUT    grid_StencilComputer_c_binding.f90 )
+
 if( atlas_HAVE_ATLAS_NUMERICS )
     generate_fortran_bindings(FORTRAN_BINDINGS ../atlas/numerics/Nabla.h)
     generate_fortran_bindings(FORTRAN_BINDINGS ../atlas/numerics/Nabla.h)
@@ -230,6 +234,7 @@ ecbuild_add_library( TARGET atlas_f
         grid/atlas_GridDistribution_module.F90
         grid/atlas_Vertical_module.F90
         grid/atlas_Partitioner_module.F90
+        grid/atlas_StencilComputer_module.F90
         mesh/atlas_MeshBuilder_module.F90
         mesh/atlas_MeshGenerator_module.F90
         mesh/atlas_Mesh_module.F90
diff --git a/src/atlas_f/atlas_module.F90 b/src/atlas_f/atlas_module.F90
index 97ddb45ae..24de358e7 100644
--- a/src/atlas_f/atlas_module.F90
+++ b/src/atlas_f/atlas_module.F90
@@ -89,6 +89,11 @@ module atlas_module
     & atlas_RegionalGrid
 use atlas_Vertical_module, only :&
     & atlas_Vertical
+use atlas_StencilComputer_module, only: &
+    & atlas_StructuredGrid_ComputeNorth, &
+    & atlas_StructuredGrid_ComputeWest, &
+    & atlas_StructuredGrid_ComputeStencil, &
+    & atlas_StructuredGrid_Stencil
 use atlas_functionspace_EdgeColumns_module, only: &
     & atlas_functionspace_EdgeColumns
 use atlas_functionspace_CellColumns_module, only: &
diff --git a/src/atlas_f/grid/atlas_StencilComputer_module.F90 b/src/atlas_f/grid/atlas_StencilComputer_module.F90
new file mode 100644
index 000000000..8ec0a1b0f
--- /dev/null
+++ b/src/atlas_f/grid/atlas_StencilComputer_module.F90
@@ -0,0 +1,339 @@
+! (C) Copyright 2013 ECMWF.
+!
+! 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.
+! In applying this licence, ECMWF does not waive the privileges and immunities
+! granted to it by virtue of its status as an intergovernmental organisation nor
+! does it submit to any jurisdiction.
+
+#include "atlas/atlas_f.h"
+
+module atlas_StencilComputer_module
+
+use fckit_shared_object_module, only: fckit_shared_object, fckit_c_deleter
+use atlas_kinds_module, only : ATLAS_KIND_IDX
+use, intrinsic :: iso_c_binding, only : c_ptr, c_int
+
+implicit none
+
+private :: fckit_shared_object, fckit_c_deleter
+private :: c_ptr, c_int
+private :: ATLAS_KIND_IDX
+
+public :: atlas_StructuredGrid_ComputeNorth
+public :: atlas_StructuredGrid_ComputeWest
+public :: atlas_StructuredGrid_ComputeStencil
+public :: atlas_StructuredGrid_Stencil
+
+private
+
+!------------------------------------------------------------------------------
+TYPE, extends(fckit_shared_object) :: atlas_StructuredGrid_ComputeNorth
+
+! Purpose :
+! -------
+!   *atlas_StructuredGrid_ComputeNorth* : To compute latitude index north of latitude
+
+! Methods :
+! -------
+
+! Author :
+! ------
+!   9-Jan-2024 Willem Deconinck     *ECMWF*
+
+!------------------------------------------------------------------------------
+contains
+
+  procedure, private :: atlas_StructuredGrid_ComputeNorth__execute_real32
+  procedure, private :: atlas_StructuredGrid_ComputeNorth__execute_real64
+  generic :: execute => atlas_StructuredGrid_ComputeNorth__execute_real32, &
+                      & atlas_StructuredGrid_ComputeNorth__execute_real64
+
+#if FCKIT_FINAL_NOT_INHERITING
+  final :: atlas_StructuredGrid_ComputeNorth__final_auto
+#endif
+
+END TYPE atlas_StructuredGrid_ComputeNorth
+
+interface atlas_StructuredGrid_ComputeNorth
+  module procedure atlas_StructuredGrid_ComputeNorth__ctor
+end interface
+
+!------------------------------------------------------------------------------
+
+!------------------------------------------------------------------------------
+TYPE, extends(fckit_shared_object) :: atlas_StructuredGrid_ComputeWest
+
+! Purpose :
+! -------
+!   *atlas_StructuredGrid_ComputeWest* : To compute longitude index west of longitude at given latitude index
+
+! Methods :
+! -------
+
+! Author :
+! ------
+!   9-Jan-2024 Willem Deconinck     *ECMWF*
+
+!------------------------------------------------------------------------------
+contains
+
+  procedure, private :: atlas_StructuredGrid_ComputeWest__execute_real32
+  procedure, private :: atlas_StructuredGrid_ComputeWest__execute_real64
+  generic :: execute => atlas_StructuredGrid_ComputeWest__execute_real32, &
+                      & atlas_StructuredGrid_ComputeWest__execute_real64
+
+#if FCKIT_FINAL_NOT_INHERITING
+  final :: atlas_StructuredGrid_ComputeWest__final_auto
+#endif
+
+END TYPE atlas_StructuredGrid_ComputeWest
+
+interface atlas_StructuredGrid_ComputeWest
+  module procedure atlas_StructuredGrid_ComputeWest__ctor
+end interface
+
+!------------------------------------------------------------------------------
+
+integer(c_int), parameter, private :: STENCIL_MAX_WIDTH = 6 !-> maximum 6x6 stencil
+
+TYPE atlas_StructuredGrid_Stencil
+  integer(ATLAS_KIND_IDX) :: width
+  integer(ATLAS_KIND_IDX) :: j_begin
+  integer(ATLAS_KIND_IDX) :: i_begin(STENCIL_MAX_WIDTH) ! on stack
+contains
+  procedure, pass :: write => atlas_StructuredGrid_Stencil__write
+  generic, public :: write(FORMATTED) => write
+  procedure, public :: i => atlas_StructuredGrid_Stencil__i
+  procedure, public :: j => atlas_StructuredGrid_Stencil__j
+END TYPE atlas_StructuredGrid_Stencil
+
+TYPE :: atlas_StructuredGrid_ComputeStencil
+  integer(c_int) :: halo
+  integer(ATLAS_KIND_IDX) :: stencil_width
+  integer(ATLAS_KIND_IDX) :: stencil_offset
+  type(atlas_StructuredGrid_ComputeNorth) :: compute_north
+  type(atlas_StructuredGrid_ComputeWest)  :: compute_west
+contains
+  procedure :: setup => atlas_StructuredGrid_ComputeStencil__setup
+  procedure :: execute => atlas_StructuredGrid_ComputeStencil__execute_real64
+
+  procedure :: assignment_operator => atlas_StructuredGrid_ComputeStencil__assignment
+  generic, public :: assignment(=) => assignment_operator
+
+  procedure :: final => atlas_StructuredGrid_ComputeStencil__final
+END TYPE
+
+! Better not use, use setup member function instead !
+!interface atlas_StructuredGrid_ComputeStencil
+!  module procedure atlas_StructuredGrid_ComputeStencil__ctor
+!end interface
+
+
+!========================================================
+contains
+!========================================================
+
+
+! -----------------------------------------------------------------------------
+! Destructor
+
+#if FCKIT_FINAL_NOT_INHERITING
+ATLAS_FINAL subroutine atlas_StructuredGrid_ComputeNorth__final_auto(this)
+  type(atlas_StructuredGrid_ComputeNorth), intent(inout) :: this
+#if FCKIT_FINAL_NOT_PROPAGATING
+  call this%final()
+#endif
+  FCKIT_SUPPRESS_UNUSED( this )
+end subroutine
+#endif
+
+#if FCKIT_FINAL_NOT_INHERITING
+ATLAS_FINAL subroutine atlas_StructuredGrid_ComputeWest__final_auto(this)
+  type(atlas_StructuredGrid_ComputeWest), intent(inout) :: this
+#if FCKIT_FINAL_NOT_PROPAGATING
+  call this%final()
+#endif
+  FCKIT_SUPPRESS_UNUSED( this )
+end subroutine
+#endif
+
+
+! -----------------------------------------------------------------------------
+! Constructors
+
+function atlas_StructuredGrid_ComputeNorth__ctor(grid, halo) result(this)
+  use, intrinsic :: iso_c_binding, only : c_int
+  use fckit_c_interop_module, only: c_str
+  use atlas_grid_StencilComputer_c_binding, only : atlas__grid__ComputeNorth__new, atlas__grid__ComputeNorth__delete
+  use atlas_grid_module, only : atlas_StructuredGrid
+  implicit none
+  type(atlas_StructuredGrid_ComputeNorth) :: this
+  type(atlas_StructuredGrid), intent(in) :: grid
+  integer(c_int), intent(in) :: halo
+  call this%reset_c_ptr( atlas__grid__ComputeNorth__new(grid%CPTR_PGIBUG_B, halo), &
+    & fckit_c_deleter(atlas__grid__ComputeNorth__delete) )
+  call this%return()
+end function
+
+function atlas_StructuredGrid_ComputeWest__ctor(grid, halo) result(this)
+  use, intrinsic :: iso_c_binding, only : c_int
+  use fckit_c_interop_module, only: c_str
+  use atlas_grid_StencilComputer_c_binding, only : atlas__grid__ComputeWest__new, atlas__grid__ComputeWest__delete
+  use atlas_grid_module, only : atlas_StructuredGrid
+  implicit none
+  type(atlas_StructuredGrid_ComputeWest) :: this
+  type(atlas_StructuredGrid), intent(in) :: grid
+  integer(c_int), intent(in) :: halo
+  call this%reset_c_ptr( atlas__grid__ComputeWest__new(grid%CPTR_PGIBUG_B, halo), &
+    & fckit_c_deleter(atlas__grid__ComputeWest__delete) )
+  call this%return()
+end function
+
+subroutine atlas_StructuredGrid_ComputeStencil__setup(this, grid, stencil_width)
+  use, intrinsic :: iso_c_binding, only : c_double
+  use atlas_grid_module, only : atlas_StructuredGrid
+  implicit none
+  class(atlas_StructuredGrid_ComputeStencil) :: this
+  type(atlas_StructuredGrid), intent(in) :: grid
+  integer(ATLAS_KIND_IDX), intent(in) :: stencil_width
+  this%stencil_width = stencil_width
+  this%halo = (stencil_width + 1) / 2
+  this%stencil_offset = stencil_width - floor(real(stencil_width,c_double) / 2._c_double + 1._c_double, ATLAS_KIND_IDX)
+  this%compute_north = atlas_StructuredGrid_ComputeNorth(grid, this%halo)
+  this%compute_west  = atlas_StructuredGrid_ComputeWest(grid, this%halo)
+end subroutine
+
+
+
+! ----------------------------------------------------------------------------------------
+
+function atlas_StructuredGrid_ComputeNorth__execute_real32(this, y) result(index)
+  use, intrinsic :: iso_c_binding, only : c_float
+  use atlas_grid_StencilComputer_c_binding, only : atlas__grid__ComputeNorth__execute_real32
+  implicit none
+  integer(ATLAS_KIND_IDX) :: index
+  class(atlas_StructuredGrid_ComputeNorth), intent(in) :: this
+  real(c_float), intent(in) :: y
+  index = atlas__grid__ComputeNorth__execute_real32(this%CPTR_PGIBUG_B, y) + 1
+end function
+
+function atlas_StructuredGrid_ComputeNorth__execute_real64(this, y) result(index)
+  use, intrinsic :: iso_c_binding, only : c_double
+  use atlas_grid_StencilComputer_c_binding, only : atlas__grid__ComputeNorth__execute_real64
+  implicit none
+  integer(ATLAS_KIND_IDX) :: index
+  class(atlas_StructuredGrid_ComputeNorth), intent(in) :: this
+  real(c_double), intent(in) :: y
+  index = atlas__grid__ComputeNorth__execute_real64(this%CPTR_PGIBUG_B, y) + 1
+end function
+! ----------------------------------------------------------------------------------------
+
+function atlas_StructuredGrid_ComputeWest__execute_real32(this, x, j) result(index)
+  use, intrinsic :: iso_c_binding, only : c_float
+  use atlas_grid_StencilComputer_c_binding, only : atlas__grid__ComputeWest__execute_real32
+  implicit none
+  integer(ATLAS_KIND_IDX) :: index
+  class(atlas_StructuredGrid_ComputeWest), intent(in) :: this
+  real(c_float), intent(in) :: x
+  integer(ATLAS_KIND_IDX), intent(in) :: j
+  index = atlas__grid__ComputeWest__execute_real32(this%CPTR_PGIBUG_B, x, j-int(1,ATLAS_KIND_IDX)) + 1
+end function
+
+function atlas_StructuredGrid_ComputeWest__execute_real64(this, x, j) result(index)
+  use, intrinsic :: iso_c_binding, only : c_double
+  use atlas_grid_StencilComputer_c_binding, only : atlas__grid__ComputeWest__execute_real64
+  implicit none
+  integer(ATLAS_KIND_IDX) :: index
+  class(atlas_StructuredGrid_ComputeWest), intent(in) :: this
+  real(c_double), intent(in) :: x
+  integer(ATLAS_KIND_IDX), intent(in) :: j
+  index = atlas__grid__ComputeWest__execute_real64(this%CPTR_PGIBUG_B, x, j-int(1,ATLAS_KIND_IDX)) + 1
+end function
+! ----------------------------------------------------------------------------------------
+
+
+subroutine atlas_StructuredGrid_ComputeStencil__execute_real64(this, x, y, stencil)
+  use, intrinsic :: iso_c_binding, only : c_double
+  implicit none
+  class(atlas_StructuredGrid_ComputeStencil), intent(in) :: this
+  real(c_double), intent(in) :: x, y
+  type(atlas_StructuredGrid_Stencil), intent(inout) :: stencil
+
+  integer(ATLAS_KIND_IDX) :: jj
+  stencil%width = this%stencil_width
+
+  stencil%j_begin = this%compute_north%execute(y) - this%stencil_offset
+  do jj = 1_ATLAS_KIND_IDX, this%stencil_width
+    stencil%i_begin(jj) = this%compute_west%execute(x, stencil%j_begin + jj - 1) - this%stencil_offset
+  enddo
+end subroutine
+! ----------------------------------------------------------------------------------------
+
+subroutine atlas_StructuredGrid_Stencil__write (stencil, unit, iotype, v_list, iostat, iomsg)
+  implicit none
+  class(atlas_StructuredGrid_Stencil), intent(in) :: stencil
+    INTEGER, INTENT(IN) :: unit
+    CHARACTER(*), INTENT(IN) :: iotype
+    INTEGER, INTENT(IN)  :: v_list(:)
+    INTEGER, INTENT(OUT) :: iostat
+    CHARACTER(*), INTENT(INOUT) :: iomsg
+    integer(ATLAS_KIND_IDX) :: jlat, jlon
+    do jlat = 1, stencil%width
+      write(unit,'(a,I0,a)',advance='no',IOSTAT=iostat) " j: ", stencil%j(jlat), "    i = "
+      do jlon = 1, stencil%width
+        write(unit,'(I3)',advance='no',IOSTAT=iostat) stencil%i(jlon,jlat)
+      enddo
+      write(0,'(a)',IOSTAT=iostat) new_line('a')
+    enddo 
+end subroutine
+
+function atlas_StructuredGrid_Stencil__j(this, j_index) result(j)
+  integer(ATLAS_KIND_IDX) :: j
+  class(atlas_StructuredGrid_Stencil), intent(in) :: this
+  integer(ATLAS_KIND_IDX) :: j_index
+  j = this%j_begin + (j_index-1)
+end function
+
+function atlas_StructuredGrid_Stencil__i(this, i_index, j_index) result(i)
+  integer(ATLAS_KIND_IDX) :: i
+  class(atlas_StructuredGrid_Stencil), intent(in) :: this
+  integer(ATLAS_KIND_IDX) :: i_index
+  integer(ATLAS_KIND_IDX) :: j_index
+  i = this%i_begin(j_index) + (i_index-1)
+end function
+
+! ----------------------------------------------------------------------------------------
+
+function atlas_StructuredGrid_ComputeStencil__ctor(grid, stencil_width) result(this)
+  use atlas_grid_module, only : atlas_StructuredGrid
+  implicit none
+  type(atlas_StructuredGrid_ComputeStencil) :: this
+  type(atlas_StructuredGrid), intent(in) :: grid
+  integer(ATLAS_KIND_IDX), intent(in) :: stencil_width
+  call this%setup(grid, stencil_width)
+end function
+
+
+subroutine atlas_StructuredGrid_ComputeStencil__assignment(this, other)
+implicit none
+class(atlas_StructuredGrid_ComputeStencil), intent(inout) :: this
+class(atlas_StructuredGrid_ComputeStencil), intent(in) :: other
+call this%final()
+write(0,*) "owners = ", other%compute_north%owners()
+this%compute_north = other%compute_north
+this%compute_west = other%compute_west
+this%stencil_width = other%stencil_width
+this%stencil_offset = other%stencil_offset
+this%halo = other%halo
+end subroutine
+
+subroutine atlas_StructuredGrid_ComputeStencil__final(this)
+  class(atlas_StructuredGrid_ComputeStencil), intent(inout) :: this
+  call this%compute_north%final()
+  call this%compute_west%final()
+end subroutine
+
+! ----------------------------------------------------------------------------------------
+
+end module atlas_StencilComputer_module