Fortran 2008 provides regular Bessel functions.

I’m looking for an implementation of spherical Bessel and Hankel functions. I need them for both first and second kind, and of arbitrary order (accepts order as input).

I have what seems to be a F66 implementation in a scan of an NRL report by Blake (1972). However, I can’t find the code elsewhere and transcribing it would be error prone.

Do you know of a library that implements these?