Skip to content

Commit

Permalink
NAG fixes. (#213)
Browse files Browse the repository at this point in the history
* NAG fixes.

Co-authored-by: Willem Deconinck <[email protected]>
  • Loading branch information
DJDavies2 and wdeconinck authored Jul 15, 2024
1 parent 8981cd1 commit b698b91
Show file tree
Hide file tree
Showing 3 changed files with 49 additions and 10 deletions.
55 changes: 47 additions & 8 deletions src/atlas_f/field/atlas_Field_module.fypp
Original file line number Diff line number Diff line change
Expand Up @@ -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

!-------------------------------------------------------------------------------
Expand Down Expand Up @@ -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"
Expand All @@ -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"
Expand All @@ -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()
Expand Down
2 changes: 1 addition & 1 deletion src/atlas_f/grid/atlas_Grid_module.F90
Original file line number Diff line number Diff line change
Expand Up @@ -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(:)
Expand Down
2 changes: 1 addition & 1 deletion src/atlas_f/mesh/atlas_MeshBuilder_module.F90
Original file line number Diff line number Diff line change
Expand Up @@ -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) )
Expand Down

0 comments on commit b698b91

Please sign in to comment.