maint/fix: move nan handling into Fortran code

CHANGELOG.md

@@ -18,6 +18,9 @@ All notable changes to this project will be documented in this file.
     and array-like inputs.
   - This should have minimal impact on users, as it is a helper function
     and behavior of the calculation routines is unchanged.
+- **FIXED** Missing values are now returned as NaN directly from the Fortran
+  wrappers instead of a small sentinel value (`9.99e-38`/`9.999e-38`) that was
+  converted to NaN in Python.
 
 ## [v0.12.0] 2025-11-26
 

pymsis/msis.py

@@ -299,10 +299,6 @@ def calculate(
             input_data[:, 7:],
         )
 
-    # The Fortran code puts 9.9e-38 in as NaN
-    # or 9.99e-38, or 9.999e-38, so lets just bound the 9s
-    output[(output >= 9.9e-38) & (output < 1e-37)] = np.nan  # noqa: PLR2004
-
     return output.reshape(*input_shape, 11)
 
 

src/wrappers/msis00.F90

@@ -12,6 +12,7 @@ subroutine pyinitswitch(switch_legacy, parmpath)
 end subroutine pyinitswitch
 
 subroutine pymsiscalc(day, utsec, lon, lat, z, sflux, sfluxavg, ap, output, n)
+    use, intrinsic :: ieee_arithmetic, only: ieee_value, ieee_quiet_nan
     implicit none
 
     integer, intent(in)        :: n
@@ -27,8 +28,12 @@ subroutine pymsiscalc(day, utsec, lon, lat, z, sflux, sfluxavg, ap, output, n)
 
     integer :: i
     real :: t(2), d(9), lon_tmp ! Temporary to swap dimensions
+    real :: nan
+
+    nan = ieee_value(1.0, ieee_quiet_nan)
 
     do i=1, n
+        ! Normalize negative longitudes into [0, 360).
         if (lon(i) < 0) then
             lon_tmp = lon(i) + 360
         else
@@ -37,12 +42,11 @@ subroutine pymsiscalc(day, utsec, lon, lat, z, sflux, sfluxavg, ap, output, n)
         call gtd7d(10000 + FLOOR(day(i)), utsec(i), z(i), lat(i), lon_tmp, &
                    utsec(i)/3600. + lon_tmp/15., sfluxavg(i), &
                    sflux(i), ap(i, :), 48, d, t)
-        ! O, H, and N are set to zero below 72.5 km
-        ! Change them to NaN instead
+        ! O, H, and N are set to zero below 72.5 km, return NaN instead
         if(z(i) < 72.5) then
-            d(2) = 9.99e-38 ! NaN
-            d(7) = 9.99e-38 ! NaN
-            d(8) = 9.99e-38 ! NaN
+            d(2) = nan
+            d(7) = nan
+            d(8) = nan
         endif
         ! These mappings are to go from MSIS00 locations to MSIS2 locations
         output(i, 1) = d(6)
@@ -54,7 +58,7 @@ subroutine pymsiscalc(day, utsec, lon, lat, z, sflux, sfluxavg, ap, output, n)
         output(i, 7) = d(5)
         output(i, 8) = d(8)
         output(i, 9) = d(9)
-        output(i, 10) = 9.99e-38 ! NaN
+        output(i, 10) = nan ! MSIS-00 does not provide NO
         output(i, 11) = t(2)
     enddo
 
@@ -77,6 +81,7 @@ end subroutine pytselec
 
 subroutine pygtd7d(day, utsec, lon, lat, z, sflux, sfluxavg, ap, output, n)
 
+use, intrinsic :: ieee_arithmetic, only: ieee_value, ieee_quiet_nan
 implicit none
 
 integer, intent(in)        :: n
@@ -92,6 +97,9 @@ subroutine pygtd7d(day, utsec, lon, lat, z, sflux, sfluxavg, ap, output, n)
 
 integer :: i
 real :: t(2), d(9), lon_tmp ! Temporary to swap dimensions
+real :: nan
+
+nan = ieee_value(1.0, ieee_quiet_nan)
 
 WRITE(6, *) 'Warning: pygtd7d is deprecated and will be removed in a future version. Use pymsiscalc instead.'
 do i=1, n
@@ -103,12 +111,11 @@ subroutine pygtd7d(day, utsec, lon, lat, z, sflux, sfluxavg, ap, output, n)
   call gtd7d(10000 + FLOOR(day(i)), utsec(i), z(i), lat(i), lon_tmp, &
               utsec(i)/3600. + lon_tmp/15., sfluxavg(i), &
               sflux(i), ap(i, :), 48, d, t)
-  ! O, H, and N are set to zero below 72.5 km
-  ! Change them to NaN instead
+  ! O, H, and N are set to zero below 72.5 km, return NaN instead
   if(z(i) < 72.5) then
-    d(2) = 9.99e-38 ! NaN
-    d(7) = 9.99e-38 ! NaN
-    d(8) = 9.99e-38 ! NaN
+    d(2) = nan
+    d(7) = nan
+    d(8) = nan
   endif
   ! These mappings are to go from MSIS00 locations to MSIS2 locations
   output(i, 1) = d(6)
@@ -120,7 +127,7 @@ subroutine pygtd7d(day, utsec, lon, lat, z, sflux, sfluxavg, ap, output, n)
   output(i, 7) = d(5)
   output(i, 8) = d(8)
   output(i, 9) = d(9)
-  output(i, 10) = 9.99e-38 ! NaN
+  output(i, 10) = nan ! MSIS-00 does not provide NO
   output(i, 11) = t(2)
 enddo
 

src/wrappers/msis2.F90

@@ -27,7 +27,8 @@ subroutine pymsiscalc(day, utsec, lon, lat, z, sflux, sfluxavg, ap, output, n)
     !       code takes the order (z, lat, lon).
     !       sflux also comes before sfluxavg
     use msis_calc, only: msiscalc
-    use msis_constants, only: rp
+    use msis_constants, only: rp, dmissing
+    use, intrinsic :: ieee_arithmetic, only: ieee_value, ieee_quiet_nan
 
     implicit none
 
@@ -43,16 +44,15 @@ subroutine pymsiscalc(day, utsec, lon, lat, z, sflux, sfluxavg, ap, output, n)
     real(kind=rp), intent(out) :: output(n, 1:11)
 
     integer :: i
-    real(kind=rp) :: lon_tmp
 
     do i=1, n
-        if (lon(i) < 0) then
-            lon_tmp = lon(i) + 360
-          else
-            lon_tmp = lon(i)
-          endif
         call msiscalc(day(i), utsec(i), z(i), lat(i), lon(i), sfluxavg(i), &
                     sflux(i), ap(i, :), output(i, 11), output(i, 1:10))
     enddo
 
+    ! The model marks missing densities with the tiny sentinel ``dmissing``.
+    ! Convert those to NaN here so callers receive an unambiguous missing value
+    ! and no downstream sentinel matching is required.
+    where (output == dmissing) output = ieee_value(1.0_rp, ieee_quiet_nan)
+
 end subroutine pymsiscalc

tests/test_regression.py

@@ -67,9 +67,9 @@ def run_input_line(line, version):
         raise ValueError("Version number incorrect")
 
     # Different units in the test file (cgs)
-    x[np.isclose(x, 9.9e-38, atol=1e-38)] = np.nan
     x[:-1] *= 1e-6
     x[5] *= 1e3
+    # Missing values come back as NaN; the reference file encodes them as 9.999e-38
     x[np.isnan(x)] = 9.999e-38
     # The output print statements are messy.
     # They technically give 4 decimal places, but only truly 3 decimal places