real y(2),yp(2),work(15),mu_atm,m_keha,a(2)
integer iwork(5),iflag
common /aaa/ g_rask,ro_atm_0,p1,p2,p3
common /bbb/ at,aa
namelist /input/pi,g_rask,r_keha,c,m_keha,p_atm, &
t_atm,mu_atm,r,nt,dt,relerr,abserr,iflag,neqn,y
external func
open(20,file="dim1_rkf45_v3.in")
open(10,file="dim1_rkf45.dat")
!*******************
read(20,input)
print input
!********************
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

yp=y
wta=0.
wt=0.
wa=0.
wg=0.
t=0.
tout=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
dx=y(1)-yp(1)
yp=y
call func(t,y,a)
wt=wt+m_keha*at*dx
wa=wa+m_keha*aa*dx
wg=wg-g_rask*dx
wta=wt+wa
etot=0.5*m_keha*y(2)**2+m_keha*g_rask*y(1)
write (10,'(10f15.6)') tout,y,etot,wg,wt,wa,wta,etot-wta

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