Enable _4 tests

CMakeLists.txt

@@ -25,10 +25,6 @@ if(BUILD_4 AND BUILD_D)
   set(kinds "4" "d")
 elseif(BUILD_4 AND NOT BUILD_D)
   set(kinds "4")
-  if(BUILD_TESTING)
-    set(BUILD_TESTING OFF)
-    message(WARNING "Disabling build testing for BUILD_4-only case")
-  endif()
 elseif(BUILD_D AND NOT BUILD_4)
   set(kinds "d")
 else()

docs/user_guide.md

@@ -33,6 +33,13 @@ in ipolates_mod. The interpolation method is chosen via the first
 argument of these routines (variable IP). Sub-options are set via the
 IPOPT array.
 
+It should be noted that some routines may behave poorly or unpredictably when
+using 4-byte reals (libip_4). For instance, there is an ATAN2 function
+used for polar stereo grids where for certain grids/coordinates, floating point
+differences between 4-byte output values (~1e-7) can be amplified into
+noticeable differences in output field values. Some applications may therefore
+benefit from the use of 8-byte reals (libip_d).
+
 ### Bilinear Interpolation
 
 Bilinear interpolation is chosen by setting IP=0.

tests/CMakeLists.txt

@@ -21,32 +21,40 @@ if(${CMAKE_Fortran_COMPILER_ID} MATCHES "^(GNU)$" AND ${CMAKE_Fortran_COMPILER_V
   set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -w -fallow-argument-mismatch -fallow-invalid-boz")
 endif()
 
-foreach(kind "d")
+foreach(kind ${kinds})
+  string(TOUPPER ${kind} kind_definition)
+
   # Test ipxwafs routines
   add_executable(test_ipxwafs_${kind} test_ipxwafs.F90)
   target_link_libraries(test_ipxwafs_${kind} PUBLIC ip::ip_${kind})
   target_link_libraries(test_ipxwafs_${kind} PUBLIC sp::sp_${kind})
+  target_compile_definitions(test_ipxwafs_${kind} PRIVATE "LSIZE=${kind_definition}")
   add_test(test_ipxwafs_${kind} test_ipxwafs_${kind})
 
   # Test earth_radius_mod.
   add_executable(test_earth_radius_${kind} test_earth_radius.F90)
   target_link_libraries(test_earth_radius_${kind} PUBLIC ip::ip_${kind})
   target_link_libraries(test_earth_radius_${kind} PUBLIC sp::sp_${kind})
+  target_compile_definitions(test_earth_radius_${kind} PRIVATE "LSIZE=${kind_definition}")
   add_test(test_earth_radius_${kind} test_earth_radius_${kind})
 
   # grib-2 tests
   add_library(test_library_grib2_${kind} input_data_mod_grib2.F90 interp_mod_grib2.F90)
   target_link_libraries(test_library_grib2_${kind} PUBLIC ip::ip_${kind})
   target_link_libraries(test_library_grib2_${kind} PUBLIC sp::sp_${kind})
+  target_compile_definitions(test_library_grib2_${kind} PRIVATE "LSIZE=${kind_definition}")
 
   add_executable(tst_gdswzd_grib2_${kind} tst_gdswzd_grib2.c)
   set_target_properties(tst_gdswzd_grib2_${kind} PROPERTIES LINKER_LANGUAGE C)
+  target_compile_definitions(tst_gdswzd_grib2_${kind} PRIVATE "LSIZE=${kind_definition}")
   add_executable(test_scalar_grib2_${kind} test_scalar_grib2.F90)
   add_executable(test_vector_grib2_${kind} test_vector_grib2.F90)
 
   target_link_libraries(tst_gdswzd_grib2_${kind} PRIVATE test_library_grib2_${kind})
   target_link_libraries(test_scalar_grib2_${kind} PRIVATE test_library_grib2_${kind})
   target_link_libraries(test_vector_grib2_${kind} PRIVATE test_library_grib2_${kind})
+  target_compile_definitions(test_scalar_grib2_${kind} PRIVATE "LSIZE=${kind_definition}")
+  target_compile_definitions(test_vector_grib2_${kind} PRIVATE "LSIZE=${kind_definition}")
 
   add_test(tst_gdswzd_c_grib2_${kind} tst_gdswzd_grib2_${kind})
 
@@ -68,7 +76,7 @@ foreach(kind "d")
   add_test(test_station_points_neighbor_budget_scalar_grib2_${kind} test_scalar_grib2_${kind} -1 6)
 
   # # vector tests
-  add_test(test_lambert_biliner_vector_grib2_${kind} test_vector_grib2_${kind} 218 0)
+  add_test(test_lambert_bilinear_vector_grib2_${kind} test_vector_grib2_${kind} 218 0)
   add_test(test_gaussian_neighbor_vector_grib2_${kind} test_vector_grib2_${kind} 127 2)
   add_test(test_latlon_bilinear_vector_grib2_${kind} test_vector_grib2_${kind} 3 0)
   add_test(test_mercator_bicubic_vector_grib2_${kind} test_vector_grib2_${kind} 8 1)
@@ -88,9 +96,11 @@ foreach(kind "d")
   add_library(test_library_grib1_${kind} input_data_mod_grib1.F90 interp_mod_grib1.F90)
   target_link_libraries(test_library_grib1_${kind} PUBLIC ip::ip_${kind})
   target_link_libraries(test_library_grib1_${kind} PUBLIC sp::sp_${kind})
+  target_compile_definitions(test_library_grib1_${kind} PRIVATE "LSIZE=${kind_definition}")
 
   add_executable(tst_gdswzd_grib1_${kind} tst_gdswzd_grib1.c)
   set_target_properties(tst_gdswzd_grib1_${kind} PROPERTIES LINKER_LANGUAGE C)
+  target_compile_definitions(tst_gdswzd_grib1_${kind} PRIVATE "LSIZE=${kind_definition}")
   add_executable(test_scalar_grib1_${kind} test_scalar_grib1.F90)
   add_executable(test_vector_grib1_${kind} test_vector_grib1.F90)
 
@@ -104,12 +114,12 @@ foreach(kind "d")
   add_test(test_gaussian_neighbor_scalar_grib1_${kind} test_scalar_grib1_${kind} 127 2)
   add_test(test_latlon_bilinear_scalar_grib1_${kind} test_scalar_grib1_${kind} 3 0)
   add_test(test_mercator_bicubic_scalar_grib1_${kind} test_scalar_grib1_${kind} 8 1)
-  add_test(test_polar_stereo_neighbor_budget_scalar_${kind} test_scalar_grib1_${kind} 212 6)
-  add_test(test_rotatedB_spectral_scalar_${kind} test_scalar_grib1_${kind} 205 4)
+  add_test(test_polar_stereo_neighbor_budget_scalar_grib1_${kind} test_scalar_grib1_${kind} 212 6)
+  add_test(test_rotatedB_spectral_scalar_grib1_${kind} test_scalar_grib1_${kind} 205 4)
   add_test(test_rotatedE_budget_scalar_grib1_${kind} test_scalar_grib1_${kind} 203 3)
 
   # vector tests
-  add_test(test_lambert_biliner_vector_grib1_${kind} test_vector_grib1_${kind} 218 0)
+  add_test(test_lambert_bilinear_vector_grib1_${kind} test_vector_grib1_${kind} 218 0)
   add_test(test_gaussian_neighbor_vector_grib1_${kind} test_vector_grib1_${kind} 127 2)
   add_test(test_latlon_bilinear_vector_grib1_${kind} test_vector_grib1_${kind} 3 0)
   add_test(test_mercator_bicubic_vector_grib1_${kind} test_vector_grib1_${kind} 8 1)

tests/input_data_mod_grib1.F90

@@ -4,6 +4,13 @@
 ! to be interpolated.
 !
 ! Kyle Gerheiser June, 2021
+
+#if (LSIZE==D)
+#define REALSIZE 8
+#elif (LSIZE==4)
+#define REALSIZE 4
+#endif
+
 module input_data_mod_grib1
   implicit none
 
@@ -25,10 +32,11 @@ module input_data_mod_grib1
 
   ! integer, parameter :: missing=b'11111111111111111111111111111111'
 
-  real, allocatable, public      :: input_data(:,:)
+  real(KIND=REALSIZE), allocatable, public      :: input_data(:,:)
+  real(KIND=REALSIZE), allocatable, public      :: input_u_data(:,:)
+  real(KIND=REALSIZE), allocatable, public      :: input_v_data(:,:)
+
   logical*1, allocatable, public :: input_bitmap(:,:)
-  real, allocatable, public      :: input_u_data(:,:)
-  real, allocatable, public      :: input_v_data(:,:)
 
 
   data input_kgds /0,  360,  180, -89500, -180000, 128,  &

tests/input_data_mod_grib2.F90

@@ -4,6 +4,13 @@
 ! to be interpolated.
 !
 ! Kyle Gerheiser June, 2021
+
+#if (LSIZE==D)
+#define REALSIZE 8
+#elif (LSIZE==4)
+#define REALSIZE 4
+#endif
+
 module input_data_mod_grib2
   implicit none
 
@@ -22,10 +29,11 @@ module input_data_mod_grib2
 
   integer, parameter :: missing=4294967296
 
-  real, allocatable, public      :: input_data(:,:)
+  real(KIND=REALSIZE), allocatable, public      :: input_data(:,:)
+  real(KIND=REALSIZE), allocatable, public      :: input_u_data(:,:)
+  real(KIND=REALSIZE), allocatable, public      :: input_v_data(:,:)
+
   logical*1, allocatable, public :: input_bitmap(:,:)
-  real, allocatable, public      :: input_u_data(:,:)
-  real, allocatable, public      :: input_v_data(:,:)
 
   data input_gdtmpl /6, 255, missing, 255, missing, 255, missing, &
        360, 180, 0, missing, -89500000, -180000000, &

tests/interp_mod_grib1.F90

@@ -1,6 +1,15 @@
 ! This is a test for the NCEPLBS-ip library.
 !
 ! Kyle Gerheiser June, 2021
+
+#if (LSIZE==D)
+#define REALSIZE 8
+#define REALSIZESTR "8"
+#elif (LSIZE==4)
+#define REALSIZE 4
+#define REALSIZESTR "4"
+#endif
+
 module interp_mod_grib1
   use ip_mod
   implicit none
@@ -53,15 +62,16 @@ subroutine interp(grid, interp_opt)
     integer                   :: ip, ipopt(20), output_kgds(200)
     integer                   :: km, ibi(1), mi, iret, i, j, ibi_scalar=0
     integer                   :: i_output, j_output, mo, no, ibo(1), ibo_scalar
-    integer                   :: num_pts_diff, which_func
+    integer                   :: num_pts_diff, which_func, ntol
 
     logical*1, allocatable    :: output_bitmap(:,:)
 
-    real, allocatable         :: output_rlat(:), output_rlon(:)
-    real, allocatable         :: output_data(:,:)
+    real(KIND=REALSIZE), allocatable         :: output_rlat(:), output_rlon(:)
+    real(KIND=REALSIZE), allocatable         :: output_data(:,:)
     real(kind=4), allocatable :: baseline_data(:,:)
     real                      :: avgdiff, maxdiff
     real(kind=4)              :: output_data4
+    real                      :: abstol
 
     integer :: grd3(200)    ! global one-degree lat/lon
     data grd3 / 0, 360, 181, 90000, 0, 128, -90000,  &
@@ -91,6 +101,14 @@ subroutine interp(grid, interp_opt)
     data grd218 /3, 614, 428, 12190, -133459, 8, -95000,  &
          12191, 12191, 0, 64, 25000, 25000, 0, 0, 0, 0, 0, 0, 255, 180*0/
 
+#if (LSIZE==D)
+    abstol=0.0001
+    ntol = 0
+#elif (LSIZE==4)
+    abstol=0.05
+    ntol = 10
+#endif
+
     select case (trim(grid))
     case ('3')
        output_kgds = grd3
@@ -189,6 +207,11 @@ subroutine interp(grid, interp_opt)
                  no, output_rlat, output_rlon, ibo_scalar, output_bitmap, output_data, iret)
         endif
 
+! Uncomment to generate new baseline file:
+!        open (13, file="grid"//trim(grid)//".opt"//trim(interp_opt)//".bin_"//REALSIZESTR, access="direct", recl=mo*4)
+!        write (13, rec=1) real(output_data, kind=4)
+!        close (13)
+
         if (iret /= 0) then
            print*,'- BAD STATUS FROM IPOLATES: ', iret
            stop 6
@@ -228,7 +251,7 @@ subroutine interp(grid, interp_opt)
         do j = 1, j_output
            do i = 1, i_output
               output_data4 = real(output_data(i,j),4)
-              if ( abs(output_data4 - baseline_data(i,j)) > 0.0001) then
+              if ( abs(output_data4 - baseline_data(i,j)) > abstol) then
                  avgdiff = avgdiff + abs(output_data4-baseline_data(i,j))
                  num_pts_diff = num_pts_diff + 1
                  if (abs(output_data4-baseline_data(i,j)) > abs(maxdiff))then
@@ -247,8 +270,8 @@ subroutine interp(grid, interp_opt)
         endif
         print*,'- AVG DIFFERENCE: ', avgdiff
 
-        if (num_pts_diff > 0) then
-           print *, "Expected 0 points > 0, found: ", num_pts_diff
+        if (num_pts_diff > ntol) then
+           print *, "# DIFFERING POINTS: ", num_pts_diff
            error stop
         endif
     enddo ! which_func
@@ -315,18 +338,19 @@ subroutine interp_vector(grid, interp_opt)
     integer                   :: km, ibi(1), mi, iret, i, j, which_func
     integer                   :: i_output, j_output, mo, no, ibo(1)
     integer                   :: ibi_scalar=0, ibo_scalar
-    integer                   :: num_upts_diff, num_vpts_diff
+    integer                   :: num_upts_diff, num_vpts_diff, ntol
 
     logical*1, allocatable    :: output_bitmap(:,:)
 
-    real, allocatable         :: output_rlat(:), output_rlon(:)
-    real, allocatable         :: output_crot(:), output_srot(:)
-    real, allocatable         :: output_u_data(:,:), output_v_data(:,:)
+    real(KIND=REALSIZE), allocatable         :: output_rlat(:), output_rlon(:)
+    real(KIND=REALSIZE), allocatable         :: output_crot(:), output_srot(:)
+    real(KIND=REALSIZE), allocatable         :: output_u_data(:,:), output_v_data(:,:)
     real                      :: avg_u_diff, avg_v_diff
     real                      :: max_u_diff, max_v_diff
     real(kind=4)              :: output_data4
     real(kind=4), allocatable :: baseline_u_data(:,:)
     real(kind=4), allocatable :: baseline_v_data(:,:)
+    real                      :: abstol
 
     integer :: grd3(200)    ! global one-degree lat/lon
     data grd3 / 0, 360, 181, 90000, 0, 128, -90000,  &
@@ -356,6 +380,14 @@ subroutine interp_vector(grid, interp_opt)
     data grd218 /3, 614, 428, 12190, -133459, 8, -95000,  &
          12191, 12191, 0, 64, 25000, 25000, 0, 0, 0, 0, 0, 0, 255, 180*0/
 
+#if (LSIZE==D)
+    abstol=0.0001
+    ntol = 0
+#elif (LSIZE==4)
+    abstol=0.05
+    ntol = 10
+#endif
+
     select case (trim(grid))
     case ('3')
        output_kgds = grd3
@@ -458,6 +490,12 @@ subroutine interp_vector(grid, interp_opt)
                  ibo_scalar, output_bitmap, output_u_data, output_v_data, iret)
         endif
 
+! Uncomment to generate new baseline file:
+!        open (13, file="grid"//trim(grid)//".opt"//trim(interp_opt)//".bin_"//REALSIZESTR, access="direct", recl=mo*4)
+!        write (13, rec=1) real(output_u_data, kind=4)
+!        write (13, rec=2) real(output_v_data, kind=4)
+!        close (13)
+
         if (iret /= 0) then
            print*,'- BAD STATUS FROM IPOLATES: ', iret
            stop 6
@@ -519,15 +557,15 @@ subroutine interp_vector(grid, interp_opt)
         do j = 1, j_output
            do i = 1, i_output
               output_data4 = real(output_u_data(i,j),4)
-              if (abs(output_data4 - baseline_u_data(i,j)) > 0.0001) then
+              if (abs(output_data4 - baseline_u_data(i,j)) > abstol) then
                  avg_u_diff = avg_u_diff + abs(output_data4-baseline_u_data(i,j))
                  num_upts_diff = num_upts_diff + 1
                  if (abs(output_data4-baseline_u_data(i,j)) > abs(max_u_diff))then
                     max_u_diff = output_data4-baseline_u_data(i,j)
                  endif
               endif
               output_data4 = real(output_v_data(i,j),4)
-              if (abs(output_data4 - baseline_v_data(i,j)) > 0.0001) then
+              if (abs(output_data4 - baseline_v_data(i,j)) > abstol) then
                  avg_v_diff = avg_v_diff + abs(output_data4-baseline_v_data(i,j))
                  num_vpts_diff = num_vpts_diff + 1
                  if (abs(output_data4-baseline_v_data(i,j)) > abs(max_v_diff))then
@@ -555,8 +593,8 @@ subroutine interp_vector(grid, interp_opt)
         endif
         print*,'- AVG DIFFERENCE: ', avg_v_diff
 
-        if (num_vpts_diff + num_upts_diff > 0) then
-           print *, "Expected 0 points > 0, found: ", num_upts_diff + num_vpts_diff
+        if (num_vpts_diff + num_upts_diff > ntol) then
+           print *, "# DIFFERING POINTS: ", num_upts_diff + num_vpts_diff
            error stop
         endif
     enddo ! which_func