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,dt
real r_keha,y(2),t,tout,relerr,abserr,work(15)
integer i,n,neqn,iflag,iwork(5)
common p1,p2,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
p1=-c_keha*s_keha/(2.*m_keha)
p2=g_rask*v_keha/m_keha

dt=1.
n=3000
t=0.
tout=0.
neqn=2
y(1)=0.
y(2)=0.
iflag=1
relerr=1.e-8
abserr=1.e-8

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
print '(3f15.5)',tout,y
stop
end

subroutine func(t,y,dy)
use COESA_module
implicit none
real*8 ro
real t,y(2),dy(2),p1,p2,g_rask,v
common p1,p2,g_rask
ro=COESA_density(dble(y(1)))
dy(1)=y(2)
v=abs(y(2))
dy(2)=-g_rask+p1*ro*v*y(2)+p2*ro
return
end