implicit none
real mu_atm,m_keha,pi,p_atm,t_atm,g_rask,c_keha,r_atm
real s_keha,v_keha,p1,p2,p3,ro_0,dt
real r_keha,y(2),t,tout,relerr,abserr,work(15)
integer i,n,neqn,iflag,iwork(5)
common p1,p2,p3,ro_0,g_rask
external func
open(10,file="dim1_rkf45.dat")
pi=3.1415926536
m_keha=15.
r_keha=1.8
mu_atm=29.
p_atm=100400.
t_atm=280.
g_rask=9.814
c_keha=0.47
r_atm=8314.

s_keha=pi*r_keha**2
v_keha=4.*pi*r_keha**3/3
ro_0=p_atm*mu_atm/(r_atm*t_atm)
p1=-c_keha*s_keha/(2.*m_keha)
p2=g_rask*v_keha/m_keha
p3=mu_atm*g_rask/(r_atm*t_atm)

dt=0.1
n=50000
t=0.
tout=0.
neqn=2
y(1)=0.
y(2)=0.
iflag=1
relerr=1.e-7
abserr=1.e-7

do i=1,n
 tout=t+dt
 call rkf45(func,neqn,y,t,tout,relerr,abserr,iflag,work,iwork)
  if (iflag /= 2) then
   iflag=2
  endif
  if(y(1)<0)then
   print *,"ERROR:  x=",y(1),"<0"
   stop
  endif
 write (10,'(15f20.5)')t,y
enddo

stop
end

subroutine func(t,y,dy)
implicit none
real t,y(2),dy(2),p1,p2,p3,ro_0,g_rask,ro,v
common p1,p2,p3,ro_0,g_rask
ro=ro_0*exp(-p3*y(1))
dy(1)=y(2)
v=abs(y(2))
dy(2)=-g_rask+p1*ro*v*y(2)+p2*ro
return
end