Skip to content

NAG fixes needed for GEOSldas #1496

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
mathomp4 merged 5 commits into from
Aug 28, 2026
Merged

NAG fixes needed for GEOSldas #1496

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
e43bf0b
nag for GEOSldas
Aug 10, 2026
27c436b
Merge branch 'develop' into feature/wjiang/nag_GEOSldas
gmao-rreichle Aug 25, 2026
ae5a390
Merge branch 'develop' into feature/wjiang/nag_GEOSldas
gmao-rreichle Aug 26, 2026
ac82b0e
Merge branch 'develop' into feature/wjiang/nag_GEOSldas
sdrabenh Aug 26, 2026
679a9dd
Merge branch 'develop' into feature/wjiang/nag_GEOSldas
sdrabenh Aug 28, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 4 additions & 1 deletion ...s_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/topography/utils_topo/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,10 @@
ecbuild_add_executable (TARGET generate_scrip_cube_topo.x SOURCES generate_scrip_cube.F90 geompack.F90)
target_link_libraries (generate_scrip_cube_topo.x PRIVATE MPI::MPI_Fortran esmf)
# CMake has an OpenMP issue with NAG Fortran: https://gitlab.kitware.com/cmake/cmake/-/issues/21280
if (NOT CMAKE_Fortran_COMPILER_ID MATCHES "NAG")
if (CMAKE_Fortran_COMPILER_ID MATCHES "NAG")
# This standalone legacy source uses byte-sized kind selectors (4 and 8).
set_source_files_properties (geompack.F90 PROPERTIES COMPILE_OPTIONS "-kind=byte")
else ()
target_link_libraries(generate_scrip_cube_topo.x PRIVATE OpenMP::OpenMP_Fortran)
endif ()

Expand Down
80 changes: 40 additions & 40 deletions ...p/GEOSsurface_GridComp/Utils/Raster/preproc/topography/utils_topo/generate_scrip_cube.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -496,11 +496,11 @@ program ESMF_GenerateCSGridDescription
! Correct for longitude periodicity
lon_w = minval(node_xy(1,:))
lon_e = maxval(node_xy(1,:))
if (abs(lon_e - lon_w) > 1.5_8*pi) then
if (abs(lon_e - lon_w) > 1.5_REAL64*pi) then
if (tmp_center_lons(i,j) < pi) then
where (node_xy(1,:) > pi) node_xy_tmp(1,:) = node_xy(1,:) - 2._8*pi
where (node_xy(1,:) > pi) node_xy_tmp(1,:) = node_xy(1,:) - 2._REAL64*pi
else
where (node_xy(1,:) < pi) node_xy_tmp(1,:) = node_xy(1,:) + 2._8*pi
where (node_xy(1,:) < pi) node_xy_tmp(1,:) = node_xy(1,:) + 2._REAL64*pi
endif
endif

Expand Down Expand Up @@ -568,8 +568,8 @@ program ESMF_GenerateCSGridDescription

if (found_degenerate) then
! build a tiny fallback polygon around the center
tiny_dlon = 1.0d-2 * pi/180._8
tiny_dlat = 1.0d-2 * pi/180._8
tiny_dlon = 1.0d-2 * pi/180._REAL64
tiny_dlat = 1.0d-2 * pi/180._REAL64
clon = tmp_center_lons(i,j) ! radians
clat = tmp_center_lats(i,j) ! radians
p1 = [clon - tiny_dlon, clat - tiny_dlat]
Expand All @@ -591,8 +591,8 @@ program ESMF_GenerateCSGridDescription
end if

! write fallback into SCRIP arrays
SCRIP_CornerLon(:,n) = modulo([p1(1),p2(1),p3(1),p4(1)]*(180._8/pi),360.0_8)
SCRIP_CornerLat(:,n) = [p1(2),p2(2),p3(2),p4(2)]*(180._8/pi)
SCRIP_CornerLon(:,n) = modulo([p1(1),p2(1),p3(1),p4(1)]*(180._REAL64/pi),360.0_REAL64)
SCRIP_CornerLat(:,n) = [p1(2),p2(2),p3(2),p4(2)]*(180._REAL64/pi)
fallback_mask(n) = .true.
failed_cells = failed_cells + 1

Expand Down Expand Up @@ -649,8 +649,8 @@ program ESMF_GenerateCSGridDescription
if (bad_corner) then
clon = tmp_center_lons(i,j)
clat = tmp_center_lats(i,j)
tiny_dlon = 1.0d-4 * pi/180._8
tiny_dlat = 1.0d-4 * pi/180._8
tiny_dlon = 1.0d-4 * pi/180._REAL64
tiny_dlat = 1.0d-4 * pi/180._REAL64
node_xy_tmp(:,1) = [modulo(clon-tiny_dlon,2*pi), clat-tiny_dlat]
node_xy_tmp(:,2) = [modulo(clon+tiny_dlon,2*pi), clat-tiny_dlat]
node_xy_tmp(:,3) = [modulo(clon+tiny_dlon,2*pi), clat+tiny_dlat]
Expand All @@ -675,8 +675,8 @@ program ESMF_GenerateCSGridDescription

clon = tmp_center_lons(i,j)
clat = tmp_center_lats(i,j)
tiny_dlon = 1.0d-4 * pi/180._8
tiny_dlat = 1.0d-4 * pi/180._8
tiny_dlon = 1.0d-4 * pi/180._REAL64
tiny_dlat = 1.0d-4 * pi/180._REAL64

p1 = [modulo(clon-tiny_dlon,2*pi), clat-tiny_dlat]
p2 = [modulo(clon+tiny_dlon,2*pi), clat-tiny_dlat]
Expand All @@ -694,8 +694,8 @@ program ESMF_GenerateCSGridDescription
endif

!Write CCW corners and POSITIVE area
SCRIP_CornerLon(:,n) = modulo([p1(1),p2(1),p3(1),p4(1)]*(180._8/pi), 360.0_8)
SCRIP_CornerLat(:,n) = [p1(2),p2(2),p3(2),p4(2)]*(180._8/pi)
SCRIP_CornerLon(:,n) = modulo([p1(1),p2(1),p3(1),p4(1)]*(180._REAL64/pi), 360.0_REAL64)
SCRIP_CornerLat(:,n) = [p1(2),p2(2),p3(2),p4(2)]*(180._REAL64/pi)

SCRIP_Area(n) = get_area_spherical_polygon(p1, p2, p3, p4)

Expand All @@ -720,8 +720,8 @@ program ESMF_GenerateCSGridDescription
swap_p = p2; p2 = p4; p4 = swap_p ! make CW
end if

SCRIP_CornerLon(:,n) = modulo([p1(1),p2(1),p3(1),p4(1)]*(180._8/pi), 360.0_8)
SCRIP_CornerLat(:,n) = [p1(2),p2(2),p3(2),p4(2)]*(180._8/pi)
SCRIP_CornerLon(:,n) = modulo([p1(1),p2(1),p3(1),p4(1)]*(180._REAL64/pi), 360.0_REAL64)
SCRIP_CornerLat(:,n) = [p1(2),p2(2),p3(2),p4(2)]*(180._REAL64/pi)

SCRIP_Area(n) = sph_tri_area_rad(p1,p2,p3) + sph_tri_area_rad(p1,p3,p4) ! steradians
if (SCRIP_Area(n) <= 0.d0) SCRIP_Area(n) = 1.d-12
Expand Down Expand Up @@ -938,7 +938,7 @@ program ESMF_GenerateCSGridDescription
end do

pair_local(1) = best_dist
pair_local(2) = real(max(0,best_idx), 8)
pair_local(2) = real(max(0,best_idx), REAL64)
call MPI_Allreduce(pair_local, pair_global, 1, MPI_2DOUBLE_PRECISION, MPI_MINLOC, mpiC, mpi_err)
global_idx_max_rrfac = int(pair_global(2))

Expand Down Expand Up @@ -1705,12 +1705,12 @@ subroutine remove_duplicate_points(points, num_points, eps, unique_points, num_u

subroutine safe_reorder_hull(node_xy_tmp, hull, p1, p2, p3, p4, n, i, j)
implicit none
real(8), intent(in) :: node_xy_tmp(2,4)
real(REAL64), intent(in) :: node_xy_tmp(2,4)
integer, intent(in) :: hull(4)
real(8), intent(out) :: p1(2), p2(2), p3(2), p4(2)
real(REAL64), intent(out) :: p1(2), p2(2), p3(2), p4(2)
integer, intent(in) :: n, i, j
integer :: idx
real(8) :: tmp_nodes(2,4)
real(REAL64) :: tmp_nodes(2,4)

do idx = 1,4
if (hull(idx) < 1 .or. hull(idx) > 4) then
Expand Down Expand Up @@ -1765,17 +1765,17 @@ subroutine sort_by_angles(n, angles, idx)

subroutine reorder_hull_quad(node_xy, p1, p2, p3, p4)
implicit none
real(8), intent(in) :: node_xy(2,4)
real(8), intent(out) :: p1(2), p2(2), p3(2), p4(2)
real(8) :: centroid(2), angles(4)
real(REAL64), intent(in) :: node_xy(2,4)
real(REAL64), intent(out) :: p1(2), p2(2), p3(2), p4(2)
real(REAL64) :: centroid(2), angles(4)
integer :: i, order(4), temp_order
real(8) :: temp_angle, temp_x(4), temp_y(4)
real(REAL64) :: temp_angle, temp_x(4), temp_y(4)
logical :: swapped
real(8) :: xyz1(3), xyz2(3), xyz3(3), xyz4(3), normal(3), xyz_centroid(3), orientation, tmp_p(2)
real(REAL64) :: xyz1(3), xyz2(3), xyz3(3), xyz4(3), normal(3), xyz_centroid(3), orientation, tmp_p(2)
logical, parameter :: verbose = .false.

! Compute centroid
centroid = [sum(node_xy(1,:))/4.0_8, sum(node_xy(2,:))/4.0_8]
centroid = [sum(node_xy(1,:))/4.0_REAL64, sum(node_xy(2,:))/4.0_REAL64]

! Compute angles from centroid to each corner
do i = 1,4
Expand Down Expand Up @@ -1829,10 +1829,10 @@ subroutine reorder_hull_quad(node_xy, p1, p2, p3, p4)
xyz4 = lonlat_to_xyz(p4(1), p4(2))

normal = cross(xyz1, xyz2) + cross(xyz2, xyz3) + cross(xyz3, xyz4) + cross(xyz4, xyz1)
xyz_centroid = (xyz1 + xyz2 + xyz3 + xyz4) / 4.0_8
xyz_centroid = (xyz1 + xyz2 + xyz3 + xyz4) / 4.0_REAL64
orientation = dot_product(normal, xyz_centroid)

if (orientation > 0.0_8) then
if (orientation > 0.0_REAL64) then
tmp_p = p2
p2 = p4
p4 = tmp_p
Expand All @@ -1851,10 +1851,10 @@ end subroutine reorder_hull_quad
! Calculate signed spherical polygon area using Girard's formula:
function get_signed_area_spherical_polygon(p1, p2, p3, p4) result(area_signed)
implicit none
real(8), intent(in) :: p1(2), p2(2), p3(2), p4(2)
real(8) :: area_signed
real(8) :: angles(4), excess
real(8), dimension(3) :: xyz1, xyz2, xyz3, xyz4
real(REAL64), intent(in) :: p1(2), p2(2), p3(2), p4(2)
real(REAL64) :: area_signed
real(REAL64) :: angles(4), excess
real(REAL64), dimension(3) :: xyz1, xyz2, xyz3, xyz4

xyz1 = lonlat_to_xyz(p1(1), p1(2))
xyz2 = lonlat_to_xyz(p2(1), p2(2))
Expand All @@ -1866,17 +1866,17 @@ function get_signed_area_spherical_polygon(p1, p2, p3, p4) result(area_signed)
angles(3) = vertex_angle(xyz2, xyz3, xyz4)
angles(4) = vertex_angle(xyz3, xyz4, xyz1)

excess = sum(angles) - 2.0_8*pi
excess = sum(angles) - 2.0_REAL64*pi
area_signed = excess ! negative for CW, positive for CCW
end function

! Helper function to calculate vertex angle:
function vertex_angle(xyzA, xyzB, xyzC) result(angle)
implicit none
real(8), intent(in) :: xyzA(3), xyzB(3), xyzC(3)
real(8) :: angle
real(8), dimension(3) :: AB, CB, cross_prod
real(8) :: norm_cross, dot_prod
real(REAL64), intent(in) :: xyzA(3), xyzB(3), xyzC(3)
real(REAL64) :: angle
real(REAL64), dimension(3) :: AB, CB, cross_prod
real(REAL64) :: norm_cross, dot_prod

AB = xyzA - xyzB
CB = xyzC - xyzB
Expand All @@ -1890,17 +1890,17 @@ function vertex_angle(xyzA, xyzB, xyzC) result(angle)

! Convert spherical coordinates (lon,lat) to Cartesian coordinates
function lonlat_to_xyz(lon, lat) result(xyz)
real(kind=8), intent(in) :: lon, lat
real(kind=8) :: xyz(3)
real(REAL64), intent(in) :: lon, lat
real(REAL64) :: xyz(3)

xyz(1) = cos(lat) * cos(lon)
xyz(2) = cos(lat) * sin(lon)
xyz(3) = sin(lat)
end function

function cross(a, b) result(c)
real(kind=8), intent(in) :: a(3), b(3)
real(kind=8) :: c(3)
real(REAL64), intent(in) :: a(3), b(3)
real(REAL64) :: c(3)

c(1) = a(2)*b(3) - a(3)*b(2)
c(2) = a(3)*b(1) - a(1)*b(3)
Expand Down