Rungekuta orden proyectil fortran

Se dispara un proyectil al aire con un angulo de 45 con respecto al suelo, con u(0)=v(0)=150 m por segundo, donde u, v son las velocidades horizontal y vertical respectivamente. las ecuaciones demovimiento están dadas por:
u' = -cVu
v' = -g-cVv
donde u y v son funciones del tiempo y

V=sqrt(u^2 + v^2)
c= coeficiente de arrastre 0.005
g= 9.81 m por segundo

A colntinuacion unprograma en fortran que calcula la trayectoria del proyectil por el método de runge kuta de segundo y tercer orden.

PROGRAM RK2

real :: t, x, y, k1, k2, l1, l2

data t, x, y, c, u,v , h/ 0.0, 0.0, 0.0, 0.005, 150, 150, 0.1/

print *, ' t u v x y'

open(1, file="datprobl_1a.dat", status="unknown")

DO n=1,200t = t+h
v1=sqrt(u**2+v**2)
!metodo de RK2
k1 = -c*v1*u*h
l1 = (-9.81-c*v1*v)*h

v2=sqrt((u+k1)**2+(v+l1)**2)

k2 = -c*v2*(u+k1)*h
l2 = (-9.81-c*v2*(v+l1))*hUU=u+(k1+k2)/2
VV=v+(l1+l2)/2

!Calculamos la trayectoria del proyectil:

x= x+ 0.5*(u+UU)*h

y= y+ 0.5*(v+VV)*h

u=UU
v=VV

write(1,*)t,x,yprint*,t,UU,VV,x,y

if(y .LT. 0) stop !Detiene el ciclo do cuando el proyectil "toca el suelo"

END DO

close(1)

END PROGRAM RK2

PROGRAM RK3

real :: t, x, y,k1, k2, l1, l2

data t, x, y, c, u, v , h/ 0.0, 0.0, 0.0, 0.005, 150, 150, 0.1/

print *, ' t u v x y'

open(1,file="datprobl_1a.dat", status="unknown")

DO n=1,200

t = t+h
v1=sqrt(u**2+v**2)

!Metodo de RK3
k1 = -c*v1*u*h
l1 = (-9.81-c*v1*v)*h

v2=sqrt((u+(0.5)*k1)**2+(v+(0.5)*l1)**2)

k2 =-c*v2*(u+k1)*(3./2)*h
l2 =(-9.81-c*v2*(v+l1))*(3./2)*h

v3=sqrt((u-k1+2*k2)**2+(v-l1+2*l2)**2)

k3=-c*v3*(u-k1+2*k2)*2*h
l3 = (-9.81-c*v3*(v-l1+2*l2))*2*h

UU=u+(1./6)*(k1+4*k2+k3)...