what is wrong with this code?can someone walk me through this?

Heres a fortran code to solve the equation of motion of very small
particle in a low reynolds number flow. It is not working in all
compilers, I have a working version of this code, but it doesnt work
for all inputs. Can someone tell me what is wrong with this and fix it
so that it works for a stokes number of 0.1 through 2. Stokes number
is calculated in the program. You need to work backwards by playing
with the input file and see for what values of input the stokes number
is 0.1 through 2.

my particle hardly moves. I have no way of even telling if my values
are correct. What should I do? Pls help.
my input file. it contains values of t_max, alpha, Rho, Rp, Rf, Uf in
that order exactly
example working input file values are:
600,0.05,2165,3.00E-07,5.20E-06,0.003 in that order.

source code:

IMPLICIT NONE
REAL*8 DERIV_FX,
DERIV_FY,N,RP,RF,RHO,UF,alpha,K,R0,U
REAL*8
Nu,x,y,y0,THETA,DT,T0,UT0,PI,T,THETA0,Stokes
REAL*8
UR0,b,x1,VX,VY,R,UT_K,UR_K,THETA2,MP,R2,x_i
INTEGER t_max, time
REAL*8 n_t
EXTERNAL FX, FY, DERIV_FX, DERIV_FY
COMMON Nu, RP, RF, MP, UF, alpha, K,
PI, Stokes, x_i
CALL read_parameters (t_max, alpha,
RHO, RP, RF, UF)
PI = 3.1415
Nu = 0.0000181
MP = RHO * ((4.0/3.0)*PI*RP**3.0)
Stokes = (2.0*RP) **2*RHO*UF/
(6*PI*Nu*2.0*RF)
write (6,*) 'The Stokes is equal to:'
write (6,*) Stokes
b = RF * alpha**(-1.0/2.0)
K =
-0.5*log(alpha)-0.75+alpha-0.25*alpha**2.0
n_t = 1.0
5 y0 = n_t*b
x= (-1.0)*(b**2.0 -
y0**2.0)**(1.0/2.0)
R0 = (x**2.0 + y0**2.0)**(1.0/2.0)
THETA0 = atan(y0/x)


THETA0= PI/2.0 + THETA0
DT = 0.000000001
T0 = 0.0
T =T0
THETA = THETA0
R = R0
Y= Y0
write (6,*) 'n_t=', n_t, 'y=', y
U = UF

UT0 = (U/(2.0*K))*(2.0*log(R/RF)+1+alpha-(RF**2.0/R**2.0) *(1.0-alpha/
2.0)-(3.0*alpha/2.0)*(R**2.0/RF**2.0))*sin(THETA)
UR0 = (U/(2.0*K))*(2.0*log(R/RF)-1.0+alpha+(RF**2.0/R**2.0) *(1.0-
alpha/2.0)-(alpha/2.0)*(R**2.0/RF**2.0))*cos(THETA)
VX= UT0
VY = UR0
OPEN (11, FILE='PARTICLEOUT.DAT')
WRITE (11,*) 'OUTPUT VARIABLES =
time,THETA, R, x, y, VT, UT,VR, UR'
DO 7 time = 1, t_max
x_i =x
CALL KUTTA(x,y,VX,VY,T,DT,FX,FY)
THETA2 = atan(y/x)
R2 = (x**2.0 + y**2.0)**(1.0/2.0)
UT_K = (U/(2.0*K))*(2.0*log(R2/RF)+1.0+alpha-(RF**2.0/R2**2.0)*(1.0-
alpha/2.0)-3.0*alpha/2.0*(R2**2.0/RF**2.0))*sin(THETA2)
UR_K = (U/(2.0*K))*(2.0*log(R2/RF)-1.0+alpha+(RF**2.0/R2**2.0)*(1.0-
alpha/2.0)-alpha/2.0*(R2**2.0/RF**2.0))*cos(THETA2)
WRITE(11,10)
time,THETA2,R2,x,y,VX,UT_K,VY,UR_K

10 FORMAT(I1,9E15.6)

if ((x**2+y**2) .LE. (RF+RP)**2) then
write (6,*) 'The Stokes is equal to:'
write (6,*) Stokes
stop
endif

if (n_t .LE. 8.0E-8) then
close (11)
write (6,*) 'Neglegible single fiber
efficiency under these conditions'

stop
endif

if ((x-x_i) .LE. 0.0) then
n_t = n_t - (n_t/200.0)
close (11)
goto 5
stop
endif
7 CONTINUE
CLOSE (11)
END

SUBROUTINE KUTTA(x,y,VX,VY,T,DT,FX,FY)
IMPLICIT NONE
REAL*8
D1X,D1Y,DT,VX,VY,D1U,D1V,D2X,D2Y,D2U,D2V
REAL*8
D3X,D3Y,D3U,D3V,D4X,D4Y,D4U,D4V,T,x,y
REAL*8 DERIV_FX, DERIV_FY

DERIV_FX = 0.0
DERIV_FY =0.0

D1X = DT * VX
D1Y = DT * VY

CALL FX(x,y,VX,DERIV_FX)
CALL FY(x,y,VY,DERIV_FY)

D1U = DT * DERIV_FX
D1V= DT * DERIV_FY

D2X = DT * (VX+D1U/2.0)
D2Y = DT * (VY+D1V/2.0)

CALL FX(x+D1X/2.0, y+D1Y/2.0, VX+D1U/
2.0, DERIV_FX)
CALL FY(x+D1X/2.0, y+D1Y/2.0, VY+D1V/
2.0, DERIV_FY)

D2U = DT * DERIV_FX
D2V = DT * DERIV_FY

D3X = DT * (VX+D2U/2.0)
D3Y = DT * (VY+D2V/2.0)

CALL FX(x+D2X/2.0, y+D2Y/2.0, VX+D2U/
2.0, DERIV_FX)
CALL FY(x+D2X/2.0, y+D2Y/2.0, VY+D2V/
2.0, DERIV_FY)

D3U = DT * DERIV_FX
D3V = DT * DERIV_FY

D4X = DT * (VX+D3U)
D4Y = DT * (VY+D3V)

CALL FX(x+D3X, y+D3Y, Vx+D3U,
DERIV_FX)
CALL FY(x+D3X, y+D3Y, VY+D3V,
DERIV_FY)

D4U = DT * DERIV_FX
D4V= DT * DERIV_FY

T = T + DT
X= x + (D1X + 2.0*D2X + 2.0*D3X +
D4X) / 6.0
y =y + (D1Y + 2.0*D2Y + 2.0*D3Y +
D4Y) / 6.0
VX = VX + (D1U + 2.0*D2U + 2.0*D3U +
D4U) / 6.0
VY = VY + (D1V + 2.0*D2V + 2.0*D3V +
D4V) / 6.0

RETURN
END

SUBROUTINE FX (x, y, VX, DERIV_FX)

IMPLICIT NONE
REAL*8
UT_K,Sk,K,Nu,alpha,y,RP,RF,MP,UF,x,VX
REAL*8
Stokes,DERIV_FX,PI,THETA,R,x_i,U
COMMON
Nu,RP,RF,MP,UF,alpha,K,PI,Stokes,x_i

R= (x**2.0 + y**2.0)**(1.0/2.0)
U = UF

THETA = atan(y/x)
UT_K = (U/(2.0*K))*(2.0*log(R/RF)
+1.0+alpha-(RF**2.0/R**2.0)*(1.0-alpha/2.0)-(3.0*alpha/2.0)*(R**2.0/
RF**2.0))*sin(THETA)
Sk = 6.0*PI*Nu*RP/MP

DERIV_FX = Sk*(UT_K-VX)

RETURN
END

SUBROUTINE FY(x, y, VY, DERIV_FY)

IMPLICIT NONE
REAL*8 UR_K, Sk, K, Nu, alpha, y, RP, RF, MP,
UF, THETA, VY
REAL*8 Stokes,DERIV_FY,PI,x,R,x_i,U
COMMON Nu,RP,RF,MP,UF,alpha,K,PI,Stokes,x_i

R= (x**2.0 + y**2.0)**(1.0/2.0)
U = UF


THETA = atan(y/x)
UR_K=( (U/(2.0*K))*(2.0*log(R/RF)-1.0+alpha+
(RF**2.0/R**2.0)*(1.0-alpha/2.0)-(alpha/2.0)*(R**2.0/
RF**2.0))*cos(THETA))

Sk = 6.0*PI*Nu*RP/MP


DERIV_FY = Sk* (UR_K-VY)

RETURN
END

SUBROUTINE read_parameters (t_max,alpha, RHO,
RP, RF, UF)
implicit none
integer t_max
real*8 alpha,RHO,RP,RF,UF
open (10, file='particle.par')

read(10,*) t_max

read(10,*) alpha

read(10,*) RHO

read(10,*) RP

read(10,*) RF

read(10,*) UF

close (10)

write (6,*) '*** Parameters read from file
particle.par'
write (6,*) '***'

return
end


.