fix: use global fortran lock instead of per-library

pymsis/msis.py

@@ -13,13 +13,13 @@
 
 # We need to point to the MSIS parameter file that was installed with the Python package
 _MSIS_PARAMETER_PATH = str(Path(__file__).resolve().parent) + "/"
+# A single global lock guarding all calls into the Fortran code.
+# per-library isn't sufficient because the Fortran code uses global state
+_lock = threading.Lock()
 for lib in [msis00f, msis20f, msis21f]:
     # Store the previous options to avoid reinitializing the model
     # each iteration unless necessary
     lib._last_used_options = None
-    # Anytime we call into the Fortran code, we need to lock
-    # to avoid threading issues
-    lib._lock = threading.Lock()
 
 
 class Variable(IntEnum):
@@ -282,7 +282,7 @@ def calculate(
                 "one of the valid version numbers: (0, 2.0, 2.1)"
             )
 
-    with msis_lib._lock:
+    with _lock:
         # Only reinitialize the model if the options have changed
         if msis_lib._last_used_options != options:
             msis_lib.pyinitswitch(options, parmpath=_MSIS_PARAMETER_PATH)

src/wrappers/msis2.F90

@@ -45,6 +45,11 @@ subroutine pymsiscalc(day, utsec, lon, lat, z, sflux, sfluxavg, ap, output, n)
 
     integer :: i
 
+    ! Zero the output before calling into the model. There are some issues with
+    ! potential nan-leakage on unitiailzed arrays in free-threading builds,
+    ! so just set it to 0 to avoid any potential issues.
+    output = 0.0_rp
+
     do i=1, n
         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))

tests/test_msis.py

@@ -508,12 +508,20 @@ def test_multithreaded(
         [f107a] * n,
         ap * n,
     )
+    expected_output20 = np.squeeze(pymsis.calculate(*input_data, version=2.0))
+    expected_output20_with_options = np.squeeze(
+        pymsis.calculate(*input_data, version=2.0, options=[0] * 25)
+    )
+
     # Create a tuple of items (version, options, expected_output)
-    # 3 items cycled over 100 times
+    # cycled over 100 times. Mixing versions 2.0 and 2.1 with differing
+    # options forces frequent re-initialization of multiple libraries.
     list_of_inputs = [
         (0, None, np.tile(expected_output00, (n, 1))),
         (2.1, None, np.tile(expected_output, (n, 1))),
         (2.1, [0] * 25, np.tile(expected_output_with_options, (n, 1))),
+        (2.0, None, expected_output20),
+        (2.0, [0] * 25, expected_output20_with_options),
     ] * 100
 
     def run_function(input_items):