real y(2),yp(2),work(15),mu_atm,m_keha,a(2)
integer iwork(5)
common /aaa/ g_rask,ro_atm_0,p1,p2,p3
common /bbb/ at,aa
external func
open(10,file="dim1_rkf45.dat")
!*******************
pi=3.1415926536
g_rask=9.814
r_keha=1.8
m_keha=15.
c=0.47
p_atm=100000.
t_atm=293.
mu_atm=29.
r=8314.
nt=10000
dt=0.01
relerr=1.e-7
abserr=1.e-7
iflag=1
t=0.
neqn=2
!********************
s_keha=pi*r_keha**2
v_keha=4.*pi*r_keha**3/3
ro_atm_0=p_atm*mu_atm/(r*t_atm)
p1=g_rask*mu_atm/(r*t_atm)
p2=0.5*c*s_keha/m_keha
p3=g_rask*v_keha/m_keha

y(1)=0.
y(2)=0.

do it=1,nt
tout=t+dt
call rkf45(func,neqn,y,t,tout,relerr,abserr,iflag,work,iwork)
iflag=2
if(y(1).le.0.) then
stop
else
call func(t,y,a)
etot=0.5*m_keha*y(2)**2+m_keha*g_rask*y(1)
write (10,'(10f15.6)') tout,y,at,aa,at+aa-g_rask,etot

endif
enddo

stop
end

subroutine func(t,y,yp)
real y(2),yp(2),m_keha
common /aaa/ g_rask,ro_atm_0,p1,p2,p3
common /bbb/ at,aa
yp(1)=y(2)
ro=ro_atm_0*exp(-p1*y(1))
v=abs(y(2))
at=-p2*ro*v*y(2)
aa=p3*ro
yp(2)=-g_rask+at+aa
return
end