      PROGRAM free_fall                     ! no air resistance
      IMPLICIT NONE
      REAL*8 a, dt, g, t, v, y
      INTEGER counter, nshow
      CALL initial(y,v,a,g,t,dt)            ! initial conditions and parameters
      CALL print_table(y,v,a,t,nshow)       ! print initial conditions
      counter = 0
      DO WHILE(y .GE. 0)
         CALL Euler(y,v,a,g,t,dt)
         counter = counter + 1
         IF(mod(counter,nshow) .EQ. 0) THEN
            CALL print_table(y,v,a,t,nshow)
         END IF
      END DO
      CALL print_table(y,v,a,t,nshow)       ! print values at surface
      END
C
      SUBROUTINE initial(y,v,a,g,t,dt)
      IMPLICIT NONE
      REAL*8 y, v, a, g, t, dt
      t = 0                                 ! initial time (sec)
      y = 10                                ! initial height (m)
      v = 0                                 ! initial velocity
      g = 9.8D0                             ! (magnitude) of accel due to gravity
      a = -g
      WRITE(*,'(A,$)') 'time step dt = '
      READ(*,*) dt
      END
C
      SUBROUTINE Euler(y,v,a,g,t,dt)
      IMPLICIT NONE
      REAL*8 y, v, a, g, t, dt
      y = y + v*dt                          ! use velocity at beginning of interval
*     following included to remind us that acceleration is constant
      a = -g                                ! y positive upward
      v = v + a*dt
      t = t + dt
      END
C
      SUBROUTINE print_table(y,v,a,t,nshow)
      IMPLICIT NONE
      REAL*8 y, v, a, t
      INTEGER nshow
      IF(t .EQ. 0) THEN
         WRITE(*,'(A,$)') 'number of time steps between output = '
         READ(*,*) nshow
         WRITE(*,*)
         WRITE(*,'(A)') 
     +     '     time (s)        y (m) velocity (m/s) accel (m/s^2)'
         WRITE(*,*)
      END IF
      WRITE(*,'(4F13.6)') t, y, v, a
      END
