diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/topography/utils_topo/CMakeLists.txt b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/topography/utils_topo/CMakeLists.txt index 828e2ed7d..f7bd45668 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/topography/utils_topo/CMakeLists.txt +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/topography/utils_topo/CMakeLists.txt @@ -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 () diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/topography/utils_topo/generate_scrip_cube.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/topography/utils_topo/generate_scrip_cube.F90 index 329270410..7035358a0 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/topography/utils_topo/generate_scrip_cube.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/topography/utils_topo/generate_scrip_cube.F90 @@ -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 @@ -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] @@ -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 @@ -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] @@ -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] @@ -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) @@ -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 @@ -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)) @@ -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 @@ -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 @@ -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 @@ -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)) @@ -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 @@ -1890,8 +1890,8 @@ 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) @@ -1899,8 +1899,8 @@ function lonlat_to_xyz(lon, lat) result(xyz) 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)