c cfall2.f c c
      module constants
         integer, parameter :: np=2000, dbl=selected_real_kind(14,100)
         real(dbl) :: g=9.807,dtmin=.001, kspg=10., mass = 100.0
c
c    np - number of elements available in solution arrays
c    dbl - KIND for real variables
c    g  -  constant of gravitational acceleration
c    dtmin -  minimum allowed time step size
c    kspg  -  spring constant for bungy cord
c    mass -  Mass of the fool jumping off of the bridge
c
      end module constants
c
      program fall
      use constants
      implicit none
c
c   Program to calculate the dynamics of a falling body
c
c   John Mahaffy    4/15/95
c
c   Arrays:
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
c   In this program, I am using allocatate just to save space in the
c   executable program file (a.out).  No attempt is made to estimate a size.
c   Module "constants" communicates between subroutines.
c
      real(dbl), allocatable :: v(:),z(:),t(:), zreal(:)
      real(dbl) dt
      integer nsteps
c
      allocate (v(np),z(np),t(np),zreal(np))
      call input(z,dt)
      call odesolve(v,z,t,dt,nsteps)
      call realans(t,z,nsteps,zreal)
      call output (t,z,zreal,v,nsteps)
      stop
      end
c
      subroutine input (z,dt)
      use constants
      implicit none
c
c   Obtain user input for initial height and time step
c
c   John Mahaffy    4/15/95
c
c  Output Arguments:
c     z(1)   -  initial height
c     dt     -  integration time step
c
      real(dbl) z(*),dt
c
      write(6,'(a)',advance='no') ' Initial height (meters): '
      read *, z(1)
      write(6,'(a)',advance='no') 'Time step size (seconds): '
      read *, dt
      if(dt.le.0.) dt=dtmin
      return
      end
c
      subroutine odesolve(v,z,t,dt,nsteps)
      use constants
      implicit none
c
c   Solve the Ordinary Differential Equation of motion for the body
c
c    John Mahaffy   5/16/95
c
c   Arguments:
c   Input
c     dt     -   timestep size
c   Output:
c     v    -   velocity
c     z    -   height
c     t    -   time
c     nsteps - last step in the integration
c
      real(dbl) v(*),z(*),t(*),dt, c0, c1
      integer i,nsteps
c
c     g  -  gravitational acceleration
c     dtmin - time step selected if 0.0 is entered
c     kspg - spring constant
c     mass - mass of falling body
c
c
c   Solve the equation for a falling body on a spring
c
c     d v              k                      d z
c     ---    =  - g + --- ( z0 - z ) ,         ---   =  v
c     d t              m                      d t
c
c
c   Set remaining initial conditions:
c
      t(1)=0.
      v(1)=0.
c
c   Set some useful constants
c
      c1=dt*kspg/mass
      c0= c1*z(1)-g*dt
c
c     Now loop through time steps until z goes negative or we run out of space
c     This time integration is FIRST ORDER accurate (error proportional to dt)
c
      do 100 i=2,np
         v(i)= v(i-1)+c0-c1*z(i-1)
         z(i)= z(i-1)+dt*.5*(v(i)+v(i-1))
         t(i)=t(i-1)+dt
         if(z(i).lt.0.) go to 200
 100     continue
      write(6,*) 'Ran out of space to continue integration'
      write(6,*) ' Last height was ',z(np),' meters'
      i=np
 200  nsteps=i
c     return
      end
c
      subroutine realans(t,z,nsteps,zreal)
      use constants
      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     John Mahaffy 5/16/95
c
c     g  -  gravitational acceleration
c     dtmin - time step selected if 0.0 is entered
c     kspg - spring constant
c     mass - mass of falling body
c
      real(dbl) t(*),z(*),zreal(*),c1,c2
      integer i,nsteps
      c1=g*mass/kspg
c     
      c2=sqrt(kspg/mass)
c        
      do 10 i=1,nsteps
c     
  10  zreal(i)=c1*cos(c2*t(i)) + z(1)-c1
c        
      return
      end
c
      subroutine output(t,z,zreal,v,nsteps)
      use constants, only : dbl
      implicit none
c
c   Outputs the full results of the time integration
c
c   John Mahaffy    4/15/95
c
      real(dbl) v(*),z(*),t(*), zreal(*)
      integer nsteps,i
      print *, 'Results are in fall.output'
      open (11,file='fall.output')
      write(11,2000)
      do 300 i=1,nsteps
      write(11,2001) t(i),z(i),zreal(i),z(i)-zreal(i)
 300  continue
 2000 format (21x,'Computed',8x,'Actual',/,
     &        6x,'Time',12x,'Height',9x,'Height',9x,'Error')
 2001 format (1x,1p,4e15.7)
      return
      end
c
c c