module realkind
           integer, parameter :: dbl=selected_real_kind(14,100)
      end module realkind
      module constants
        use realkind
        real(dbl) :: g=9.807_dbl,dtmin=.005_dbl,m=10.0_dbl,cd=.40_dbl,
     &               rho=1.15_dbl,a=.008_dbl,v0=300._dbl, trmax=50._dbl
c
c     g -  gravitational acceleration
c     dtmin  - minimum permitted time step
c     m -  mass of the shell
c     cd - drag force coefficient
c     rho - air density
c     a   - shell cross-sectional area
c     v0  - muzzle velocity
c     trmax - maximum permitted time in the solution
c
      end module constants
      module arrays
        use realkind
        integer, parameter :: lenarr=10000
        real(dbl) :: v(lenarr),z(lenarr),t(lenarr), zreal(lenarr)
c
c     v   -   velocity at each time integration step
c     z   -   height at  each time integration step
c     t   -   time for each corresponding v and z
c     zreal - Actual height at time t(i) for comparison with computed z
c     lenarr - length of each array
c
c     Note: normally v, z, t, and zreal would be allocatable and lenarr
c           would be an integer variable obtained later based on selected
c           time step and maximum permitted problem time.  However a bug
c           in our version of dbx prevents debugging in this situation.
c
c
      endmodule arrays
c
      program shell
      use realkind
      implicit none
c
c   Program to calculate the dynamics of a vertical cannon shot
c
c     dt     -  integration time step
c     nsteps   - number of steps taken in the time solution
c     tmax     - time of maximum height
c     zmax     - maximum height
c     thit     - time shell hits the ground
c
      real(dbl) dt,tmax,zmax,thit
      integer nsteps
c
      call input(dt)
      call odesolve(dt,nsteps)
      call realans(nsteps,tmax,zmax,thit)
      call printarr(nsteps)
      call timpact(nsteps,tmax,zmax,thit)
      stop
      end
c
      subroutine input (dt)
      use arrays
      use constants
      implicit none
      real(dbl) dt
c
c   Obtain user input
c
c  Output Arguments:
c
c     dt     -  integration time step
c
c
      write(*,*) ' What time step size do you want (seconds)?'
      read *, dt
      if(dt.le.dtmin) then
         write(*,*) 'Time step must be greater than', dtmin, ' s.'
         write(*,*) 'Time step reset to the minimum allowed value.'
         dt=dtmin
      endif
c
c     Set the initial conditions for the problem
c
      t(1) = 0.
      v(1) = v0
      z(1) = 0.
c
      return
      end
c

      subroutine odesolve(dt,nsteps)
      use constants
      use arrays
      implicit none
c
c   Solve the Ordinary Differential Equation of motion for the body
c   Calculates velocity and height (v and z) at each time step.
c
c   Arguments:
c   Input
c     dt     -   timestep size
c   Output:
c     nsteps - last step in the integration
c
      real(dbl) dt,rdt,c
      real (dbl) a1,b1,c1,a0
      integer nsteps,i
c
c   Solve the equation for a falling body
c
c     d v                           d z
c     ---    =  - g - C |v|v        ---   =  v
c     d t                           d t
c
c   Set the an air friction coefficient
c
      c=cd*rho*a/(2*m)
c
c     Now loop through time steps until z goes negative or we run out of space
c
      rdt=.5_dbl/dt
      do 100 i=2,lenarr
c        v(i)= v(i-1)-dt*(g + c*abs(v(i-1))*v(i-1))
         a0=dt*c*sign(1._dbl,v(i-1))
         a1=a0*.25_dbl
         b1=1._dbl+.5_dbl*a0*v(i-1)
         c1=dt*g-v(i-1) +a1*v(i-1)**2
         v(i)=(-b1+sqrt(b1**2-4*a1*c1))/(2*a1)
         z(i)= z(i-1)+dt*.5_dbl*(v(i)+v(i-1))
         t(i)=t(i-1)+dt
         if(z(i).lt.0.) go to 200
 100     continue
      write(*,*) 'Impact not obtained before last time step'
      write(*,*) ' Last height was ',z(lenarr),' meters'
 200  nsteps=i
c     return
      end
c
      subroutine realans(nsteps,tmax,zmax,thit)
      use constants
      use arrays
      implicit none
c
c     Get the values of the analytic solution to the differential equation
c     for each time point to check the numerical accuracy.
c
c     nsteps   - number of steps taken in the time solution
c     tmax     - time of maximum height
c     zmax     - maximum height
c     thit     - time shell hits the ground
c
      real(dbl) c,b,b1
      real(dbl) tmax,zmax,thit
      integer nsteps,i
c
c   Set the an air friction coefficient
c
      c=cd*rho*a/(2*m)
c
      b=sqrt(g*c)
      tmax=atan(sqrt(c/g)*v(1))/b
      zmax= -log(cos(tmax*b))/c
      do 10 i=1,nsteps
      if(v(i).gt.0.0) then
      	zreal(i)=-log(cos(tmax*b)/cos(b*(tmax-t(i))))/c
      else
        zreal(i)=zmax -log(cosh(b*(tmax-t(i))))/c
      endif
  10  continue
      b1=exp(c*zmax)
      thit=tmax+log(b1+sqrt(b1**2-1))/b
      return
      end
c
      subroutine printarr(nsteps)
      use arrays
      use constants
      implicit none
c
c   Outputs the full results of the time integration
c
      integer nsteps,i
      open (11,file='shell.output')
      write(11,2000)
      do 300 i=1,nsteps
      write(11,2001) t(i),v(i),z(i),zreal(i)
 300  continue
 2000 format (31x,'Computed',8x,'Actual',/,
     &        6x,'Time',9x,'Velocity', 8x,'Height',8x,'Height')
 2001 format (1x,1p,4e15.7)
      return
      end
c
      subroutine timpact(nsteps,tmax,zmax,thit)
c
c    Calculate and output the time of impact and velocity impact
c
      use constants
      use arrays
      implicit none
c
c     nsteps   - number of steps taken in the time solution
c     tmax     - time of maximum height
c     zmax     - maximum height
c     thit     - time shell hits the ground
c     vhit     - velocity of impact
c     vterm    - terminal velocity (must be > vhit)
c
      real(dbl) c,b,b1
      real(dbl) tmax,zmax,thit,thitcalc,vhit,vterm,wx
      integer nsteps
c
      if (z(nsteps).gt.0.0) then
         write(*,*) 'Calculation ended before impact'
         return
         endif
      wx=z(nsteps-1)/(z(nsteps-1)-z(nsteps))
      thitcalc=t(nsteps-1)+wx*(t(nsteps)-t(nsteps-1))
      vhit=v(nsteps-1)+wx*(v(nsteps)-v(nsteps-1))
      vterm=sqrt(2.*m*g/(cd*rho*a))
      write (6,*)' Calculated Impact Velocity = ', vhit
      write (6,*)' Terminal Velocity          = ', vterm
      write (6,*)' Calculated Impact Time     = ', thitcalc
      write (6,*)' Theoretical Impact Time    = ', thit
      write (6,*)' Theoretical Time of Max z  = ', tmax
      write (6,*)' Theoretical Maximum Height = ', zmax
      return
      end
c