MODULE Shared
  USE, INTRINSIC :: ISO_Fortran_env, dp=>REAL64 !modern DOUBLE PRECISION
  IMPLICIT NONE
  INTEGER :: i, i_max, i_max_read
  INTEGER, PARAMETER :: xt_dat = 10 !manage UNIT numbers as a resource
  REAL(dp) :: x_read, v_read, const0, const0_read, gamma, gamma_read, A_0, A_0_read, omega, omega_read, dt, dt_read
  REAL(dp), DIMENSION(:), ALLOCATABLE :: x, v, a, E, dE, t
  REAL(dp), PARAMETER :: m = 1.0
  CHARACTER(*), PARAMETER :: my_format = '(2f22.15)' !.15 is dp
END MODULE Shared

PROGRAM driven
  USE Shared
  IMPLICIT NONE
  EXTERNAL Verlet
  WRITE(6, '(a34)', advance='no') "Initial position of the mass [m]: "
  READ *, x_read
  WRITE(6, '(a38)', advance='no') "Initial velocity of the mass [m/sec]: "
  READ *, v_read
  WRITE(6, '(a14)', advance='no') "Value of k/m: "
  READ *, const0_read
  const0 = const0_read
  WRITE(6, '(a43)', advance='no') "Value of the damping coefficient [kg/sec]: "
  READ *, gamma_read
  gamma = gamma_read
  WRITE(6, '(a36)', advance='no') "Amplitude of the driving force [N]: "
  READ *, A_0_read
  A_0 = A_0_read
  WRITE(6, '(a37)', advance='no') "Frequency of the driving force [Hz]: "
  READ *, omega_read
  omega = omega_read
  WRITE(6, '(a31)', advance='no') "Time-step of the system [sec]: "
  READ *, dt_read
  dt = dt_read
  WRITE(6, '(a7)', advance='no') "Steps: "
  READ *, i_max_read
  i_max = i_max_read
  WRITE(6, '(a23, f4.1, a4)') "Total simulation time: ", FLOAT(i_max) * dt, " sec"
  ALLOCATE(x(i_max + 1), v(i_max + 1), a(i_max + 1), E(i_max + 1), dE(i_max + 1), t(i_max + 1))
  x(1) = x_read
  v(1) = v_read
  t(1) = 0.0_dp
  i = 1
  damped: DO
    IF (i > i_max) EXIT
    CALL Verlet
    i = i + 1
  END DO damped
  OPEN(xt_dat, file='xt.dat', access='sequential', status='unknown', action='write')
  WRITE(xt_dat, my_format) (t(i), x(i), i=1, i_max)
  CLOSE(xt_dat)
  DEALLOCATE(x, v, a, E, dE, t)
END PROGRAM driven

SUBROUTINE Verlet()
  USE Shared
  IMPLICIT NONE
  REAL(dp) :: Force, x_temp, v_temp, t_temp
  EXTERNAL Force
  x_temp = x(i)
  v_temp = v(i)
  t_temp = t(i)
  x_temp = x(i)
  v_temp = v(i)
  a(i) = Force(x_temp, v_temp, t_temp) / m
  x(i + 1) = x(i) + v(i) * dt + 0.5_dp * a(i) * dt * dt
  x_temp = x(i + 1)
  v_temp = v(i)
  t(i + 1) = t(i) + dt
  t_temp = t(i + 1)
  a(i + 1) = Force(x_temp, v_temp, t_temp) / m
  v(i + 1) = v(i) + 0.5_dp * (a(i + 1) + a(i)) * dt
  x_temp = x(i + 1)
  v_temp = v(i + 1)
END SUBROUTINE Verlet

REAL(dp) FUNCTION Force(x_temps, v_temps, t_temps)
  USE Shared
  IMPLICIT NONE
  REAL(dp) :: x_temps, v_temps, t_temps
  REAL(dp) :: f_t
  f_t = A_0 * COS(omega * t_temps)
  Force = -const0 * x_temps - gamma * v_temps + f_t
END FUNCTION Force