implicit none
real y(2),work(15),mu,m_keha,pi,r_keha,s_keha,c_keha
real v_keha,p_atm,t_atm,r_atm,m_maa,r_maa,g_const,p1,p2,p3
real relerr,abserr,t,tout,dt,rho_atm_0,c,ro_atm_0,g_rask
integer iwork(5),nt,i,neqn,iflag

common g_const,r_maa,m_maa,p1,p2,p3,rho_atm_0
external func
open(10,file="rkf45.dat")

pi=3.1415926536
r_keha=0.01
m_keha=15.
c_keha=0.47
g_const=6.6743015e-11
m_maa=5.97219e+24
mu=29.
r_atm=8314.
t_atm=293.
p_atm=100000.
r_maa=6371000.
s_keha=pi*r_keha**2
v_keha=4.*pi*r_keha**3/3
rho_atm_0=p_atm*mu/(r_atm*t_atm)

dt=0.01
nt=400000
y(1)=0.
y(2)=10.
call rask(g_const,m_maa,r_maa,0.,g_rask)

neqn=2
t=0.
tout=0.
relerr=1.e-7
abserr=1.e-7
iflag=1
p1=0.5*c_keha*s_keha/m_keha
p2=mu*g_rask/(r_atm*t_atm)
p3=v_keha/m_keha

do i=1,nt
tout=t+dt
call rkf45(func,neqn,y,t,tout,relerr,abserr,iflag,work,iwork)
iflag=2

write (10,*)tout,y

enddo
stop
end

subroutine func(t,y,yp)
implicit none
real y(2),yp(2),p1,p2,p3,m_maa,r_maa
real rho,rho_atm_0,g_const,t,g_rask
common g_const,r_maa,m_maa,p1,p2,p3,rho_atm_0
yp(1)=y(2)
call rask(g_const,m_maa,r_maa,y(1),g_rask)
rho=rho_atm_0*exp(-p2*y(1)/(1.+y(1)/r_maa))
yp(2)=-g_rask-p1*rho*abs(y(2))*y(2)+g_rask*rho*p3
return
end

subroutine rask(g_const,m_maa,r_maa,h,g_rask)
implicit none
real g_const,m_maa,r_maa,g_rask,h
g_rask=g_const*m_maa/(r_maa+h)**2
return
end