implicit none
real y(4),work(27),mu_atm,m_keha,m_maa,a(4)
real t,tout,relerr,abserr,v0,alpha,c_keha,r_keha,p_atm
real r_gaas,s_keha,ro_atm_0,gm,t_atm,p1,p2,p3,pi,h,h0
real v_keha,g_rask,dt,g_const,aa,ro,r_maa,vv,r
integer iwork(5),i,nt,neqn,iflag
common gm,g_const,r_maa,m_maa,p2,p3
external func

open(20,file="maandumine.dat")
pi=3.1415926536
c_keha=0.47
r_keha=2.
m_keha=2800.
g_const=6.67408e-11
m_maa=5.9722e+24
gm=g_const*m_maa*1.e-9
r_maa=6371.

s_keha=pi*r_keha**2
v_keha=4.*pi*r_keha**3/3
p2=s_keha*c_keha*500./m_keha
p3=v_keha*1.e-3/m_keha

neqn=4
dt=1.
nt=200000
abserr=1.e-9
relerr=1.e-9
v0=7.7884079
alpha=0.75
h0=200.
      
y(1)=r_maa+h0
y(2)=0.
y(3)=-v0*sin(pi*alpha/180.)
y(4)=v0*cos(pi*alpha/180.)
t=0
tout=0
iflag=1
    
do i=1,nt
tout=t+dt
call rkf45(func,neqn,y,t,tout,relerr,abserr,iflag,work,iwork)

if(iflag.ne.2) then
iflag=2
endif

vv=sqrt(y(3)**2+y(4)**2)
r=sqrt(y(1)**2+y(2)**2)
h=r-r_maa
call func(t,y,a)
aa=sqrt(a(3)**2+a(4)**2)

if (h.le.5.)then
r_keha=4.
s_keha=pi*r_keha**2
v_keha=4.*pi*r_keha**3/3.
p2=c_keha*s_keha*500./m_keha
p3=v_keha*1.e-3/m_keha
endif

if(h.le.0.)then
stop
endif

write(20, '(10f15.5)') tout,y,h,vv*1000,aa*1000

enddo

stop
end

subroutine func(t,y,dy)
use COESA_module
implicit none
!use iso_fortran_env, wp => real64
!real(wp) :: hc,ro
real*8 hc,ro
real y(4),dy(4),gm,r_maa,p2,p3
real h,rr,vv,g_rask,m_maa,t,g_const
common gm,g_const,r_maa,m_maa,p2,p3
rr=sqrt(y(1)**2+y(2)**2)
vv=sqrt(y(3)**2+y(4)**2)
h=rr-r_maa

!****** COESA subroutine for air density calculation ******
hc=h*1000.
ro=COESA_density(hc)
!**********************************************************

call rask(g_const,r_maa,m_maa,h,g_rask)
dy(1)=y(3)
dy(2)=y(4)
dy(3)=-y(1)*gm/rr**3-p2*ro*vv*y(3)+p3*g_rask*ro*y(1)/rr
dy(4)=-y(2)*gm/rr**3-p2*ro*vv*y(4)+p3*g_rask*ro*y(2)/rr
return
end

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