      PROGRAM sho                           ! simple harmonic oscillator
      IMPLICIT NONE
      REAL*8 dt, dt2, t, v, w2, x
      INTEGER counter, nshow, kbhit
      LOGICAL Ltmp1_
      CALL initial(x,v,w2,t,dt,dt2,nshow)
      counter = 0
      Ltmp1_ = .TRUE.
      DO WHILE(Ltmp1_)
         CALL update(x,v,w2,t,dt,dt2)
         counter = counter + 1
         IF(counter .EQ. nshow) THEN
            CALL output(x,v,t)
            counter = 0
         END IF
         Ltmp1_ = kbhit() .EQ. 0
      END DO
      END
C
      SUBROUTINE initial(x0,v,w2,t,dt,dt2,nshow)
      IMPLICIT NONE
      REAL*8 show_time, x0, v, w2, t, dt, dt2
      INTEGER nshow
      WRITE(*,'(A,$)') 'initial position = '
      READ(*,*) x0                          ! meters
      t = 0
      v = 0                                 ! initial velocity
*     natural (angular) frequency
      WRITE(*,'(A,$)') 'ratio of k/m = '
      READ(*,*) w2
      WRITE(*,'(A,$)') 'time step = '
      READ(*,*) dt
      dt2 = 0.5D0*dt
      show_time = 0.1D0                     ! time interval between output
      nshow = int(show_time/dt)
      WRITE(*,'(A)') '         time     position     velocity'
      WRITE(*,*)                            ! skip line
      END
C
      SUBROUTINE update(x,v,w2,t,dt,dt2)
*     Euler-Richardson algorithm
      IMPLICIT NONE
      REAL*8 a, amid, vmid, xmid, x, v, w2, t, dt, dt2
      a = -w2*x
      vmid = v + a*dt2                      ! velocity at midpoint
      xmid = x + v*dt2                      ! position at midpoint
      amid = -w2*xmid                       ! acceleration at midpoint
      v = v + amid*dt
      x = x + vmid*dt
      t = t + dt
      END
C
      SUBROUTINE output(x,v,t)
*     - print number with leading space or minus sign
*     % print leading zeroes as '0'
      IMPLICIT NONE
      REAL*8 x, v, t
      WRITE(*, '(3F13.4)') t, x, v
      END
