      PROGRAM free_fall
*     no air resistance
      IMPLICIT NONE
      DOUBLE PRECISION y,v,a,g,t,dt
      INTEGER counter,nshow
*     initial conditions and parameters
      CALL initial(y,v,a,g,t,dt)
*     print initial conditions   
      CALL print_table(y,v,a,t,nshow)
      counter = 0 
      DO WHILE (y .GE. 0.0D0)
         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
*     print values at surface
      CALL  print_table(y,v,a,t,nshow)
      END

      SUBROUTINE initial(y,v,a,g,t,dt)
      IMPLICIT NONE
      DOUBLE PRECISION y,v,a,g,t,dt
*     initial time (sec)
      t = 0.0D0
*     initial height (m)
      y = 10.0D0
      v = 0.0D0
*     (magnitude) of accel due to gravity
      g = 9.8D0
      a = -g
      WRITE(*,*)'time step dt = '
      READ(*,*) dt
      END

      SUBROUTINE Euler(y,v,a,g,t,dt)
      IMPLICIT NONE
      DOUBLE PRECISION y,v,a,g,t,dt
*     use velocity at beginning of interval
      y = y + v*dt
*     following included to remind us that acceleration constant
*     y positive upward
      a = -g
      v = v + a*dt
      t = t + dt
      END

      SUBROUTINE print_table(y,v,a,t,nshow)
      IMPLICIT NONE
      DOUBLE PRECISION y,v,a,t
      INTEGER nshow
      IF (t .EQ. 0.0D0) THEN
         WRITE(*,*) 'number of time steps between output = '
         READ(*,*) nshow
         WRITE(*,*)
         WRITE(*,10)
         WRITE(*,*)
      END IF
      WRITE(*,20) t,y,v,a
10    FORMAT(5x,'time',7x,'y',6x,'velocity',6x,'acceleration')
20    FORMAT(3f10.3,5x,f10.3)
      END