      PROGRAM free_fall                     ! 空気抵抗なし
      IMPLICIT NONE
      REAL*8 a, dt, g, t, v, y
      INTEGER counter, nshow
      CALL initial(y,v,a,g,t,dt)            ! 初期条件とパラメータの設定
      CALL print_table(y,v,a,t,nshow)       ! 初期条件を出力する
      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)       ! 地面での値を出力する
      END
C
      SUBROUTINE initial(y,v,a,g,t,dt)
      IMPLICIT NONE
      REAL*8 y, v, a, g, t, dt
      t = 0                                 ! 時刻の初期値 (sec)
      y = 10                                ! 高さの初期値 (m)
      v = 0                                 ! 初速度
      g = 9.8D0                             ! 重力加速度(の大きさ)
      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                          ! 区間の初めでの速度を使う
*     加速度が一定であることを忘れないように次の行を入れておく
      a = -g                                ! y 軸は上方が正
      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
