C
      PROGRAM drag
*     assume drag force proportional to v^2
      IMPLICIT NONE
      REAL*8 dt, dt_2, g, t, v, vt2, y, y0
      INTEGER counter, nshow
      CALL initial(y,y0,v,t,g,vt2,dt,dt_2)
      CALL print_table(y,y0,t,dt,nshow)
      DO WHILE(.NOT.(t .GE. 0))
         CALL Euler_Richardson(y,v,t,g,vt2,dt,dt_2)
      END DO
      CALL print_table(y,y0,t,dt,nshow)     ! values at t = 0
      counter = 0
      DO WHILE(t .LE. 0.80D0)
         CALL Euler_Richardson(y,v,t,g,vt2,dt,dt_2)
         counter = counter + 1
         IF(mod(counter,nshow) .EQ. 0) THEN
            CALL print_table(y,y0,t,dt,nshow)
         END IF
      END DO
      END
C
      SUBROUTINE initial(y,y0,v,t0,g,vt2,dt,dt_2)
      IMPLICIT NONE
      REAL*8 vt, y, y0, v, t0, g, vt2, dt, dt_2
      t0 = -0.132D0                         ! initial time (see Table 3.1)
      y0 = 0
      y = y0
      v = 0                                 ! initial velocity
      WRITE(*,'(A,$)') 'terminal velocity (m/s) = '
      READ(*,*) vt
      vt2 = vt*vt
      dt = 0.001D0
      dt_2 = 0.5D0*dt
      g = 9.8D0
      END
C
      SUBROUTINE Euler_Richardson(y,v,t,g,vt2,dt,dt_2)
      IMPLICIT NONE
      REAL*8 a, amid, v2, vmid, vmid2, ymid, y, v, t, g, vt2, dt, dt_2
      v2 = v*v
      a = -g*(1 - v2/vt2)
      vmid = v + a*dt_2                     ! velocity at midpoint
      ymid = y + v*dt_2                     ! position at midpoint
      vmid2 = vmid*vmid
      amid = -g*(1 - vmid2/vt2)             ! acceleration at midpoint
      v = v + amid*dt
      y = y + vmid*dt
      t = t + dt
      END
C
      SUBROUTINE print_table(y,y0,t,dt,nshow)
      IMPLICIT NONE
      REAL*8 distance_fallen, show_time, y, y0, t, dt
      INTEGER nshow
      IF(t .LT. 0) THEN
         show_time = 0.1D0                  ! time interval between output
*        choice of dt = 0.001 convenient
         nshow = int(show_time/dt)
         WRITE(*,*)
         WRITE(*,'(A)') '         time  displacement'
         WRITE(*,*)
      END IF
      distance_fallen = y0 - y
      WRITE(*,'(2F13.6)') t, distance_fallen
      END