diff --git a/src/atlas_f/field/atlas_Field_module.fypp b/src/atlas_f/field/atlas_Field_module.fypp
index c841a9a00..acaadc3a3 100644
--- a/src/atlas_f/field/atlas_Field_module.fypp
+++ b/src/atlas_f/field/atlas_Field_module.fypp
@@ -183,10 +183,17 @@ subroutine array_c_to_f_${dtype}$_r${rank}$(array_cptr,rank,shape_cptr,strides_c
   enddo
   eshape(rank) = shape(rank)
   call c_f_pointer ( array_cptr , tmp , shape=eshape )
-  #{if rank == 1}# array_fptr => tmp(1:shape(1)) #{endif}#
-  #{if rank == 2}# array_fptr => tmp(1:shape(1),1:shape(2)) #{endif}#
-  #{if rank == 3}# array_fptr => tmp(1:shape(1),1:shape(2),1:shape(3)) #{endif}#
-  #{if rank == 4}# array_fptr => tmp(1:shape(1),1:shape(2),1:shape(3),1:shape(4)) #{endif}#
+  if (associated(tmp)) then
+    #{if rank == 1}# array_fptr => tmp(1:shape(1)) #{endif}#
+    #{if rank == 2}# array_fptr => tmp(1:shape(1),1:shape(2)) #{endif}#
+    #{if rank == 3}# array_fptr => tmp(1:shape(1),1:shape(2),1:shape(3)) #{endif}#
+    #{if rank == 4}# array_fptr => tmp(1:shape(1),1:shape(2),1:shape(3),1:shape(4)) #{endif}#
+  else
+    #{if rank == 1}# allocate(array_fptr(0)) #{endif}#
+    #{if rank == 2}# allocate(array_fptr(0,0)) #{endif}#
+    #{if rank == 3}# allocate(array_fptr(0,0,0)) #{endif}#
+    #{if rank == 4}# allocate(array_fptr(0,0,0,0)) #{endif}#
+  endif
 end subroutine
 
 !-------------------------------------------------------------------------------
@@ -403,7 +410,7 @@ function atlas_Field__wrap_name_${dtype}$_r${rank}$(name,data) result(field)
   use fckit_c_interop_module, only : c_str
   type(atlas_Field) :: field
   character(len=*), intent(in) :: name
-  ${ftype}$, intent(in) :: data(${dim[rank]}$)
+  ${ftype}$, intent(in), target :: data(${dim[rank]}$)
   integer(c_int) :: shapef(${rank}$)
   integer(c_int) :: stridesf(${rank}$)
 #:if ftype != "logical"
@@ -420,12 +427,13 @@ function atlas_Field__wrap_name_${dtype}$_r${rank}$(name,data) result(field)
   call field%return()
 end function
 function atlas_Field__wrap_${dtype}$_r${rank}$(data) result(field)
+  use, intrinsic :: iso_c_binding
+  use :: fckit_c_interop_module
   use atlas_field_c_binding
   use, intrinsic :: iso_c_binding, only : c_int, c_long, c_float, c_double
-  use fckit_c_interop_module, only : c_str
   use fckit_array_module, only : array_strides, array_view1d
   type(atlas_Field) :: field
-  ${ftype}$, intent(in) :: data(${dim[rank]}$)
+  ${ftype}$, intent(in), target :: data(${dim[rank]}$)
   integer(c_int) :: shapef(${rank}$)
   integer(c_int) :: stridesf(${rank}$)
 #:if ftype != "logical"
@@ -436,7 +444,38 @@ function atlas_Field__wrap_${dtype}$_r${rank}$(data) result(field)
     data1d => array_view1d( data, int(0,c_int) )
 #:endif
   shapef = shape(data)
-  stridesf = array_strides(data)
+  if( size(data)>0 ) then
+  !!! See issue https://github.com/ecmwf/atlas/pull/213 : problem with array_strides(data) with NAG compiler.
+#:if rank == 4
+    stridesf(1) = &
+      int(c_ptr_to_loc(c_loc(data(2,1,1,1)))-c_ptr_to_loc(c_loc(data(1,1,1,1))),c_int32_t)/int(8,c_int32_t)
+    stridesf(2) = &
+      int(c_ptr_to_loc(c_loc(data(1,2,1,1)))-c_ptr_to_loc(c_loc(data(1,1,1,1))),c_int32_t)/int(8,c_int32_t)
+    stridesf(3) = &
+      int(c_ptr_to_loc(c_loc(data(1,1,2,1)))-c_ptr_to_loc(c_loc(data(1,1,1,1))),c_int32_t)/int(8,c_int32_t)
+    stridesf(4) = &
+      int(c_ptr_to_loc(c_loc(data(1,1,1,2)))-c_ptr_to_loc(c_loc(data(1,1,1,1))),c_int32_t)/int(8,c_int32_t)
+#:elif rank == 3
+    stridesf(1) = &
+      int(c_ptr_to_loc(c_loc(data(2,1,1)))-c_ptr_to_loc(c_loc(data(1,1,1))),c_int32_t)/int(8,c_int32_t)
+    stridesf(2) = &
+      int(c_ptr_to_loc(c_loc(data(1,2,1)))-c_ptr_to_loc(c_loc(data(1,1,1))),c_int32_t)/int(8,c_int32_t)
+    stridesf(3) = &
+      int(c_ptr_to_loc(c_loc(data(1,1,2)))-c_ptr_to_loc(c_loc(data(1,1,1))),c_int32_t)/int(8,c_int32_t)
+#:elif rank == 2
+    stridesf(1) = &
+      int(c_ptr_to_loc(c_loc(data(2,1)))-c_ptr_to_loc(c_loc(data(1,1))),c_int32_t)/int(8,c_int32_t)
+    stridesf(2) = &
+      int(c_ptr_to_loc(c_loc(data(1,2)))-c_ptr_to_loc(c_loc(data(1,1))),c_int32_t)/int(8,c_int32_t)
+#:elif rank == 1
+    stridesf(1) = &
+      int(c_ptr_to_loc(c_loc(data(2)))-c_ptr_to_loc(c_loc(data(1))),c_int32_t)/int(8,c_int32_t)
+#:else
+    stridesf = array_strides(data)
+#:endif
+  else
+    stridesf = 0
+  endif
   field = atlas_Field__cptr( &
     atlas__Field__wrap_${ctype}$_specf( c_str(""),data1d,size(shapef),shapef, stridesf) )
 call field%return()
diff --git a/src/atlas_f/grid/atlas_Grid_module.F90 b/src/atlas_f/grid/atlas_Grid_module.F90
index 42dd9e458..39168f74d 100644
--- a/src/atlas_f/grid/atlas_Grid_module.F90
+++ b/src/atlas_f/grid/atlas_Grid_module.F90
@@ -554,7 +554,7 @@ function atlas_UnstructuredGrid__ctor_points( xy ) result(this)
   use fckit_array_module, only : array_strides, array_view1d
   use atlas_grid_unstructured_c_binding
   type(atlas_UnstructuredGrid) :: this
-  real(c_double), intent(in) :: xy(:,:)
+  real(c_double), intent(in), target :: xy(:,:)
   integer(c_int) :: shapef(2)
   integer(c_int) :: stridesf(2)
   real(c_double), pointer :: xy1d(:)
diff --git a/src/atlas_f/mesh/atlas_MeshBuilder_module.F90 b/src/atlas_f/mesh/atlas_MeshBuilder_module.F90
index 5b1fce4dd..7d0f6f4a5 100644
--- a/src/atlas_f/mesh/atlas_MeshBuilder_module.F90
+++ b/src/atlas_f/mesh/atlas_MeshBuilder_module.F90
@@ -87,7 +87,7 @@ function atlas_TriangularMeshBuilder__build(this, &
   real(c_double), contiguous, intent(in) :: x(:), y(:), lon(:), lat(:)
   integer, intent(in) :: nb_triags
   integer(ATLAS_KIND_GIDX), intent(in) :: triag_global_index(nb_triags)
-  integer(ATLAS_KIND_GIDX), intent(in) :: triag_nodes(3,nb_triags)
+  integer(ATLAS_KIND_GIDX), intent(in), target :: triag_nodes(3,nb_triags)
 
   integer(ATLAS_KIND_GIDX), pointer :: triag_nodes_1d(:)
   triag_nodes_1d => array_view1d( triag_nodes, int(0,ATLAS_KIND_GIDX) )