From 062e78cfa54283e1cb93ed778d120c2ae166b920 Mon Sep 17 00:00:00 2001 From: Greg Lucas Date: Sat, 27 Jun 2026 21:07:01 -0600 Subject: [PATCH] maint/fix: move nan handling into Fortran code The previous nan checks were handled within the numpy-layer, instead we can reference the msis-specific dmissing quantity from within Fortran and then pass back a nan to Python directly rather than having our fuzzy check for floating point values. --- CHANGELOG.md | 3 +++ pymsis/msis.py | 4 ---- src/wrappers/msis00.F90 | 31 +++++++++++++++++++------------ src/wrappers/msis2.F90 | 14 +++++++------- tests/test_regression.py | 2 +- 5 files changed, 30 insertions(+), 24 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 74d10d6..a6a5e5b 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/pymsis/msis.py b/pymsis/msis.py index 25896fd..3588460 100644 --- a/pymsis/msis.py +++ b/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) diff --git a/src/wrappers/msis00.F90 b/src/wrappers/msis00.F90 index dd79066..0b27a5b 100644 --- a/src/wrappers/msis00.F90 +++ b/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 diff --git a/src/wrappers/msis2.F90 b/src/wrappers/msis2.F90 index 9224d37..959bf89 100644 --- a/src/wrappers/msis2.F90 +++ b/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 diff --git a/tests/test_regression.py b/tests/test_regression.py index b76c494..af8a1f3 100644 --- a/tests/test_regression.py +++ b/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