Difference between test values in 00 Fortran outputs and PyMSIS

[alokshenoy]

When the MSIS00 fortran code is compiled using gfortran -std=legacy -ffixed-form -w NRLMSISE-00.FOR -o nrlmsise

and executed using ./nrlmsise, we get outputs that (based on the fortran code) are used to validate the model outputs.

Running the same inputs through pymsis however, produce very different results for the 00 model.

e.g.

import pymsis

dates = ["2023-06-21T08:03:20"]
data = pymsis.calculate(
dates, lons=[-70], lats=[0], alts=[400] ,f107s=150, f107as=150, aps=[4], version=0
)
data

output: 
array([[1.9556999e-12, 2.7302258e+12, 4.0951235e+10, 6.7588294e+13,
        2.4679893e+12, 1.2903808e+11, 8.6020760e+07, 6.3844516e+11,
        1.2644463e+09,           nan, 8.8434930e+02]], dtype=float32)

[greglucas]

One of the major differences is unit conversions. The Fortran code by default outputs g/cm^3 and #/cm^3, while pymsis output kg/m^3 and #/m^3. Here is a longer example script to iterate through all versions of MSIS and output the rows in the same order as that screenshot with the conversion factors applied. These values are all pretty close at a quick glance and pretty consistent across versions.

import numpy as np
import pymsis

d = np.datetime64("2023") + np.timedelta64(172, "D") + np.timedelta64(75000, "s")
dates = [d]
for version in (0, 2, 2.1):
    print(f"\nTesting MSIS version {version}")
    data = pymsis.calculate(
        dates, lons=[-70], lats=[60], alts=[1000], f107s=150, f107as=150, aps=[4], version=version
    )
    for variable in [pymsis.Variable.TEMPERATURE,
                     pymsis.Variable.HE,
                     pymsis.Variable.O,
                     pymsis.Variable.N2,
                     pymsis.Variable.O2,
                     pymsis.Variable.AR,
                     pymsis.Variable.H,
                     pymsis.Variable.N,
                     pymsis.Variable.ANOMALOUS_O,
                     pymsis.Variable.MASS_DENSITY]:
        if variable is pymsis.Variable.MASS_DENSITY:
            data[0][variable] *= 1e-3
        elif variable is not pymsis.Variable.TEMPERATURE:
            data[0][variable] *= 1e-6
        print(f"{variable.name}: {data[0][variable]:.3e}")

Potential difference

The test Fortran code supplies LST manually, which is the local apparent solar time. pymsis calculates this for the user's and does not provide the ability to tune that knob. In fact, we follow the advice within the Fortran code (line 65) for calculating things in a physically realistic way.

pymsis/src/wrappers/msis00.F90

Line 38 in cdd4043

utsec(i)/3600. + lon_tmp/15., sfluxavg(i), &

C     NOTES ON INPUT VARIABLES: 
C        UT, Local Time, and Longitude are used independently in the
C        model and are not of equal importance for every situation.  
C        For the most physically realistic calculation these three
C        variables should be consistent (STL=SEC/3600+GLONG/15).
C        The Equation of Time departures from the above formula
C        for apparent local time can be included if available but
C        are of minor importance.

[greglucas]

@alokshenoy were you able to investigate this more to figure out if there is an issue or not?