Has anyone in the pangeo orbit ported the ocean neutral density calculation code to Python, or have a python wrapper for the old FORTRAN program?
I have the Python version of the Matlab “gamma NP” but I would not recommend you to use that. I never really had the time to test it properly.
Eric Firing has the Fortran wrapper here: pygamma: log
oh perfect! I had totally forgotten about that. Very easy to wrap with apply_ufunc.
I once spent some time looking closely at the original Jacket and McDougal gamma_n code. (I think it was in Fortran.) I was trying to understand why it was so slow. I remember discovering that it opened a read a file with a reference T/S profile for every single profile in the inputs. So if you passed in arrays of T / S / P with shape
Nt, Nz, Ny, Nx, it would open and read the same file
Nt * Ny * Nx times. This was hard coded! An exponentially faster approach would have been to just load that data into memory one time when the library was imported.
In any case, if anyone has a fast implementation of neutral density, we would be happy to stick in fastjmd95: