more aggressive nan check in the src

Author: JanStreffing

src/gen_ic3d.F90

@@ -20,8 +20,18 @@ MODULE g_ic3d
    USE g_comm_auto
    USE g_support
    USE g_config, only: dummy, ClimateDataPath, use_cavity
+   USE, INTRINSIC :: ISO_C_BINDING, only: C_DOUBLE, C_INT
 
    IMPLICIT NONE
+   
+   ! C-binding interface for isnan - immune to Fortran compiler flags like -ffast-math
+   INTERFACE
+      PURE FUNCTION c_isnan(x) BIND(C, name="isnan")
+         IMPORT :: C_DOUBLE, C_INT
+         REAL(C_DOUBLE), VALUE, INTENT(IN) :: x
+         INTEGER(C_INT) :: c_isnan
+      END FUNCTION c_isnan
+   END INTERFACE
 
    include 'netcdf.inc'
 
@@ -65,6 +75,30 @@ MODULE g_ic3d
 
 !============== NETCDF ==========================================
 CONTAINS
+   ! Helper function: returns .true. if x is NaN, using C library (immune to -ffast-math)
+   PURE LOGICAL FUNCTION is_nan(x)
+      REAL(KIND=8), INTENT(IN) :: x
+      is_nan = (c_isnan(x) /= 0)
+   END FUNCTION is_nan
+   
+   ! Check if any of the 4 corners at any depth has NaN
+   LOGICAL FUNCTION corners_have_nan(arr, i1, i2, j1, j2, nk)
+      REAL(KIND=8), INTENT(IN) :: arr(:,:,:)
+      INTEGER, INTENT(IN) :: i1, i2, j1, j2, nk
+      INTEGER :: ii, jj, kk
+      corners_have_nan = .false.
+      do kk = 1, nk
+         do jj = j1, j2
+            do ii = i1, i2
+               if (is_nan(arr(ii,jj,kk))) then
+                  corners_have_nan = .true.
+                  return
+               end if
+            end do
+         end do
+      end do
+   END FUNCTION corners_have_nan
+   
    SUBROUTINE nc_readGrid(partit)
    ! Read time array and grid from nc file
       IMPLICIT NONE
@@ -380,11 +414,11 @@ SUBROUTINE getcoeffld(tracers, partit, mesh)
          ncdata(nc_Nlon,:,:)=ncdata(2,:,:)
 
          ! replace nan (or fillvalue) by dummy value
-         ! Use x/=x test in addition to ieee_is_nan (ieee_is_nan fails with -ffast-math)
+         ! Use C library isnan() which is immune to Fortran -ffast-math optimizations
          do k=1,nc_Ndepth
             do j=1,nc_Nlat
                do i=1,nc_Nlon
-                  if (ncdata(i,j,k) /= ncdata(i,j,k) .or. ieee_is_nan(ncdata(i,j,k)) .or. (ncdata(i,j,k)==FILL_VALUE)) then
+                  if (is_nan(ncdata(i,j,k)) .or. (ncdata(i,j,k)==FILL_VALUE)) then
                      ncdata(i,j,k) = dummy
                   elseif (ncdata(i,j,k) < -0.99_WP*dummy .or. ncdata(i,j,k) > dummy) then
                      ! and in case the input data has other conventions on missing values:
@@ -422,9 +456,10 @@ SUBROUTINE getcoeffld(tracers, partit, mesh)
          if (x<0.)   x=x+360.
          if (x>360.) x=x-360.
          if ( min(i,j)>0 ) then
-         ! Skip if any corner has dummy OR NaN (NaN comparisons are always FALSE, so check explicitly)
-         if (any(ncdata(i:ip1,j:jp1,1) > dummy*0.99_WP) .or. &
-             any(ncdata(i:ip1,j:jp1,1) /= ncdata(i:ip1,j:jp1,1))) cycle
+         ! Skip if any corner at ANY depth has dummy OR NaN
+         ! Use C library isnan() which is immune to Fortran -ffast-math optimizations
+         if (any(ncdata(i:ip1,j:jp1,:) > dummy*0.99_WP) .or. &
+             corners_have_nan(ncdata, i, ip1, j, jp1, nc_Ndepth)) cycle
             x1 = nc_lon(i)
             x2 = nc_lon(ip1)
             y1 = nc_lat(j)
@@ -438,9 +473,9 @@ SUBROUTINE getcoeffld(tracers, partit, mesh)
                     ncdata(i,jp1,:) > 0.99_WP*dummy .OR. ncdata(ip1,jp1,:) > 0.99_WP*dummy)
                 data1d(:)=dummy
             end where
-            ! Also catch any NaN using explicit loop (WHERE doesn't work with NaN due to compiler optimizations)
+            ! Also catch any NaN using C library isnan (immune to -ffast-math)
             do k=1, nc_Ndepth
-               if (data1d(k) /= data1d(k)) data1d(k) = dummy
+               if (is_nan(data1d(k))) data1d(k) = dummy
             end do   
 
             !___________________________________________________________________
@@ -557,11 +592,10 @@ SUBROUTINE do_ic3d(tracers, partit, mesh)
          end where
 
          !_________________________________________________________________________
-         ! Sanitize NaN using explicit DO loop (WHERE doesn't work with NaN due to
-         ! compiler optimizations like -ffast-math). Uses x/=x test: only NaN satisfies this.
+         ! Sanitize NaN using C library isnan (immune to -ffast-math optimizations)
          do n=1, partit%myDim_nod2d + partit%eDim_nod2D
             do i=1, mesh%nl-1
-               if (tracers%data(current_tracer)%values(i,n) /= tracers%data(current_tracer)%values(i,n)) then
+               if (is_nan(tracers%data(current_tracer)%values(i,n))) then
                   tracers%data(current_tracer)%values(i,n) = 0.0_WP
                end if
             end do

src/gen_support.F90

@@ -8,7 +8,17 @@ module g_support
   use g_comm_auto
   use o_ARRAYS
   use g_config, only: dummy
+  USE, INTRINSIC :: ISO_C_BINDING, only: C_DOUBLE, C_INT
   implicit none
+  
+  ! C-binding interface for isnan - immune to Fortran compiler flags like -ffast-math
+  INTERFACE
+     PURE FUNCTION c_isnan_support(x) BIND(C, name="isnan")
+        IMPORT :: C_DOUBLE, C_INT
+        REAL(C_DOUBLE), VALUE, INTENT(IN) :: x
+        INTEGER(C_INT) :: c_isnan_support
+     END FUNCTION c_isnan_support
+  END INTERFACE
 
   private
   public :: smooth_nod, smooth_elem, integrate_nod, integrate_elem, extrap_nod, omp_min_max_sum1, omp_min_max_sum2, point_in_polygon
@@ -419,10 +429,10 @@ subroutine extrap_nod3D(arr, partit, mesh)
     call exchange_nod(arr, partit)
 
     !___________________________________________________________________________
-    ! First convert any NaN to dummy (NaN comparisons are always FALSE, breaks logic below)
+    ! First convert any NaN to dummy using C library isnan (immune to -ffast-math)
     do nz=1, nl-1
        do n=1, myDim_nod2D+eDim_nod2D
-          if (arr(nz,n) /= arr(nz,n)) arr(nz,n) = dummy
+          if (c_isnan_support(arr(nz,n)) /= 0) arr(nz,n) = dummy
        end do
     end do
 
@@ -449,8 +459,8 @@ subroutine extrap_nod3D(arr, partit, mesh)
                 ! loop over local vertices n 
                 do n=1, myDim_nod2D!+eDim_nod2D
                     ! found node n that has to be extrapolated
-                    ! Note: NaN comparisons are always FALSE, so also check x/=x
-                    if ( ((work_array(n)>0.99_WP*dummy) .or. (work_array(n)/=work_array(n))) .and.  (nlevels_nod2D(n)>nz)) then
+                    ! Use C library isnan (immune to -ffast-math)
+                    if ( ((work_array(n)>0.99_WP*dummy) .or. (c_isnan_support(work_array(n))/=0)) .and.  (nlevels_nod2D(n)>nz)) then
                         cnt=0
                         val=0._WP