Merge branch 'hotfix/0.38.1'

* hotfix/0.38.1:
  Update Changelog
  Version 0.38.1
  NAG fixes. (#213)
  Wrap out-of-bounds source cells into the domain (#211)

CHANGELOG.md

@@ -7,6 +7,12 @@ This project adheres to [Semantic Versioning](http://semver.org/spec/v2.0.0.html
 
 ## [Unreleased]
 
+## [0.38.1] - 2024-07-15
+
+### Fixed
+- Compilation and running with NAG compiler
+- Wrap out-of-bounds source cells into the domain (#211)
+
 ## [0.38.0] - 2024-06-20
 
 ### Added
@@ -545,6 +551,7 @@ Fix StructuredInterpolation2D with retry for failed stencils
 ## 0.13.0 - 2018-02-16
 
 [Unreleased]: https://github.com/ecmwf/atlas/compare/master...develop
+[0.38.1]: https://github.com/ecmwf/atlas/compare/0.38.0...0.38.1
 [0.38.0]: https://github.com/ecmwf/atlas/compare/0.37.0...0.38.0
 [0.37.0]: https://github.com/ecmwf/atlas/compare/0.36.0...0.37.0
 [0.36.0]: https://github.com/ecmwf/atlas/compare/0.35.1...0.36.0

VERSION

@@ -1 +1 @@
-0.38.0
+0.38.1

src/atlas/interpolation/method/unstructured/UnstructuredBilinearLonLat.cc

@@ -419,15 +419,22 @@ Method::Triplets UnstructuredBilinearLonLat::projectPointToElements(size_t ip, c
             }
             // shift cells on the east-west periodic boundary from the east to the west
             // so that the quad surrounds a point with output longitude in [0,360)
-            if (lons[0] > lons[1] || lons[3] > lons[2]) {
-                lons[0] -= 360;
-                lons[3] -= 360;
+            double minlon = std::numeric_limits<double>::max();
+            for ( int i = 0; i < 4; i++ ) {
+                minlon = std::min( minlon, lons[i] );
+            }
+            for ( int i = 0; i < 4; i++ ) {
+                if ( (lons[i] - minlon) > 180 ) {
+                    lons[i] -= 360;
+                }
             }
 
             element::Quad2D quad(
                 PointLonLat{lons[0], (*ilonlat_)(idx[0], LAT)}, PointLonLat{lons[1], (*ilonlat_)(idx[1], LAT)},
                 PointLonLat{lons[2], (*ilonlat_)(idx[2], LAT)}, PointLonLat{lons[3], (*ilonlat_)(idx[3], LAT)});
 
+            ATLAS_ASSERT( quad.validate() );
+
             if (itc == elems.begin()) {
                 inv_dist_weight_quad(quad, o_loc, inv_dist_w);
             }

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()

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(:)

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) )