Fortran - Uso de funciones externas...

Hola a todos, estoy aprendiendo a utilizar el fortran 90/95 y estoy confrontando serios problemas con las funciones. Trabajo actualmente en el compilador Silver Frost FNT95 de distribución gratuita. Mi problema surge en el llamado a la función, no se si estoy haciendo las referencias correctamente o no, el programa que envio a continuación, debe crear un archivo llamado bessel.dat, el cual muestre las tablas de datos de las funciones de Bessel de primera clase (las Jn(x)), Dicho programa me compila muy bien, pero a la hora de ejecutarlo me resenta errores de ejecución: Aqui les envio el programita:

***********************************************************************************************
program Ejemplo27
implicit none

! El programa muestra el uso de las funciones de Bessel de primera especie

! Declaración de las variables.
integer :: x ! Argumento de las funciones de Bessel de primera especie.
! real :: bessj
integer :: xx

open(unit = 5, file = "besselJ.dat", status = "unknown", action = "write", &
position = "rewind")

! Entrada de los parámetros.
print *, 'Introducir el argumento maximo a trabajar (x <= 50).'
read (*,*) xx
print *

! Realización de los cálculos e impresión de los resultados.
segundo: do x = 1, xx ! Controla el argumento de la función de Bessel.
write(5,10) x, bessj0(real(x)), bessj1(real(x)), bessj(2,real(x))
!, bessj(3,real(x)), &
! bessj(4,real(x)), bessj(5,real(x)), bessj(6,real(x))
10 format(I4, 3(2X,F11.5))
end do segundo

print *, 'Archivo creado, ejecucion finalizada.'
close(5)

contains

! ***************************************
! **** Definiciones de las funciones ****
! ***************************************

real function bessj0(x) !Returns the Bessel function J0(x) for any real x.
real :: bessj0
real, intent(in) :: x
real :: ax, xx, z
real(kind(8)):: p1, p2, p3, p4, p5, q1, q2, q3, q4, q5, r1, r2, r3, r4, &
r5, r6, s1, s2, s3, s4, s5, s6, y

! Parámetros constantes pertenecientes a J0(x).
p1 = 1.0; p2 = -0.1098628627e-2; p3 = 0.2734510407e-4; p4 = -0.2073370639e-5
p5 = 0.2093887211e-6; q1 = -0.1562499995e-1; q2 = 0.1430488765e-3
q3 = -0.6911147651e-5; q4 = 0.7621095161e-6; q5 = -0.934945152e-7
r1 = 57568490574.0; r2 = -13362590354.0; r3 = 651619640.7d0
r4 = -11214424.18d0; r5 = 77392.33017d0; r6 = -184.9052456d0
s1 = 57568490411.d0; s2 = 1029532985.d0; s3 = 9494680.718d0
s4 = 59272.64853d0; s5 = 267.8532712d0; s6 = 1.d0

if(abs(x) < 8.0)then
y = x**2
bessj0 = (r1 + y * (r2 + y * (r3 + y* (r4 + y * (r5 + y * r6)))))/(s1 + &
y * (s2 + y * (s3 + y * (s4 + y * (s5 + y * s6)))))
else
ax = abs(x)
z = 8.0/ax
y = z**2
xx = ax - 0.785398164
bessj0 = sqrt(0.636619772/ax) * (cos(xx) * (p1 + y * (p2 + y * (p3 + &
y * (p4 + y * p5)))) - z * sin(xx) * (q1 + y * (q2 + y * (q3 + &
y * (q4 + y * q5)))))
endif
return
end function bessj0

! *****************************************************************************

real function bessj1(x) !Returns the Bessel function J1(x) for any real x.
real :: bessj1
real, intent(in) ::x
real :: ax, xx, z
real(kind(8)) :: p1, p2, p3, p4, p5, q1, q2, q3, q4, q5, r1, r2, r3, r4, &
r5, r6, s1, s2, s3, s4, s5, s6, y

! Parámetros constantes pertenecientes a J1(x).
r1 = 72362614232.d0; r2 = -7895059235.d0; r3 = 242396853.1d0
r4 = -2972611.439d0; r5 = 15704.48260d0; r6 = -30.16036606d0
s1 = 144725228442.d0; s2 = 2300535178.d0; s3 = 18583304.74d0
s4 = 99447.43394d0; s5 = 376.9991397d0; s6 = 1.0d0
p1 = 1.0d0; p2 = 0.183105e-2; p3 = -0.3516396496e-4; p4 = 0.2457520174e-5
p5 = -.240337019e-6; q1 = 0.04687499995d0; q2 = -0.2002690873e-3
q3 = 0.8449199096e-5; q4 = -0.88228987e-6; q5 = 0.105787412d-6

if(abs(x) < 8.0) then
y = x**2
bessj1 = x * (r1 + y * (r2 + y * (r3 + y * (r4 + y * (r5 + y * r6)))))/ &
(s1 + y * (s2 + y * (s3 + y * (s4 + y * (s5 + y * s6)))))
else
ax = abs(x)
z = 8.0/ax
y = z**2
xx = ax - 2.356194491
bessj1 = sqrt(0.636619772/ax) * (cos(xx) * (p1 + y * (p2 + y * (p3 + y * &
(p4 + y * p5)))) - z * sin(xx) * (q1 + y * (q2 + y * (q3 + y * &
(q4 + y * q5))))) * sign(1.0,x)
end if
return
end function bessj1

! *****************************************************************************

real function bessj(n,x) ! Returns the Bessel function Jn(x) for any real
! x and n >= 2.
integer :: n
real, intent(in) :: x
real :: bessj, coef1, coef2, suma = 0.0
integer :: k, p ! Contador

coef1 = (x/2.0)**n
do k = 0, 50
p = n + k
coef2 = (-x**2/4.0)**k
suma = suma + coef2/(factrl(k) * factrl(p))
end do
bessj = coef1 * suma
return
end function bessj

! *****************************************************************************

function factrl(n) !Returns the value n! as a floating-point number.
integer :: n
real :: factrl
integer :: j, ntop
real :: a(33)!, gammln
save ntop, a
data ntop, a(1)/0,1./
if (n < 0) then
pause 'negative factorial in factrl'
else if(n <= ntop) then
factrl = a(n + 1)
else if(n <= 32) then !Fill in table up to desired value.
do j = ntop + 1, n
a(j + 1) = j * a(j)
end do
ntop = n
factrl = a(n + 1)
else
factrl = exp(gammln(n + 1.0))
endif
return
end function factrl

! ****************************************************************************

function gammln(xx) !Returns the value ln[(xx)] for xx > 0.
real :: gammln
real, intent(in) :: xx
integer :: j
real(kind(8)) :: ser, stp, tmp, x, y, cof(6)

stp = 2.5066282746310005
cof(1) = 76.18009172947146; cof(2) = -86.50532032941677
cof(3) = 24.01409824083091; cof(4) = -1.231739572450155
cof(5) = 0.1208650973866179e-2; cof(6) = -0.5395239384953e-5
x = xx
y = x
tmp = x + 5.50
tmp = (x + 0.50) * log(tmp) - tmp
ser = 1.000000000190015

do j = 1, 6
y = y + 1.0
ser = ser + cof(j)/y
end do
gammln = tmp + log(stp * ser/x)
return
end function gammln
end program Ejemplo27
********************************************************
Como podran observar, estoy haciendo uso de algunas funciones contenidas en el NUMERICAL RECIPES IN FORTRAN.

Me gustaría que me ayudaran en la correcta ejecución del programa, éste funciona satisfactoriamente cuando trabajamos con las funciones de Bessel J0(x) y J1(x), pero cuando trabajamos con las funciones J2(x) en adelante, se presentan los errores. ¿Me podrían ayudar en ello por favor?

Agradeciendo la ayuda de antemano, se despide de ustedes
Alan.