MODULE Shared
  USE, INTRINSIC :: ISO_Fortran_env, dp=>REAL64 !modern DOUBLE PRECISION
  IMPLICIT NONE
  INTEGER :: i, j
  INTEGER, PARAMETER :: i_max = 100000, max_dat = 10 !manage UNIT numbers as a resource
  REAL(dp) :: k1x, k1y, k1z, k2x, k2y, k2z, k3x, k3y, k3z, k4x, k4y, k4z, sigma, sigma_read, rho, rho_read, beta, beta_read, dt,   &
    dt_read, x_read, x_temp, y_read, y_temp, z_read, z_temp
  REAL(dp), DIMENSION(:), ALLOCATABLE :: f_x, x, f_y, y, f_z, z, z_m
  CHARACTER(*), PARAMETER :: my_format = '(2f20.15)' !.15 is dp
END MODULE Shared

PROGRAM discrete
  USE Shared
  IMPLICIT NONE
  ALLOCATE(f_x(i_max + 1), x(i_max + 1), f_y(i_max + 1), y(i_max + 1), f_z(i_max + 1), z(i_max + 1), z_m(i_max + 1))
  PRINT *, " "
  WRITE(6, '(a6)', advance='no') "  x = "
  READ *, x_read
  x(1) = x_read
  IF (x(1) > 20.0_dp .OR. x(1) < 1.0_dp) STOP
  WRITE(6, '(a6)', advance='no') "  y = "
  READ *, y_read
  y(1) = y_read
  IF (y(1) > 20.0_dp .OR. y(1) < 1.0_dp) STOP
  WRITE(6, '(a6)', advance='no') "  z = "
  READ *, z_read
  z(1) = z_read
  IF (z(1) > 20.0_dp .OR. z(1) < 1.0_dp) STOP
  WRITE(6, '(a10)', advance='no') "  sigma = "
  READ *, sigma_read
  sigma = sigma_read
  IF (sigma > 10.0_dp .OR. sigma < 1.0_dp) STOP
  WRITE(6, '(a8)', advance='no') "  rho = "
  READ *, rho_read
  rho = rho_read
  IF (rho > 30.0_dp .OR. rho < 1.0_dp) STOP
  WRITE(6, '(a9)', advance='no') "  beta = "
  READ *, beta_read
  beta = beta_read
  IF (beta > 3.0_dp .OR. beta < 0.1_dp) STOP
  WRITE(6, '(a7)', advance='no') "  dt = "
  READ *, dt_read
  dt = dt_read
  IF (dt > 0.100_dp .OR. dt < 0.001_dp) STOP
  PRINT *, " "
  i = 1
  j = 2
  OPEN(max_dat, file='max.dat', access='sequential', status='unknown', action='write')
  lorenz: DO !common fourth-order Runge-Kutta method
    IF (i > i_max - 1) EXIT !do not over count
    x_temp = x(i)
    y_temp = y(i)
    z_temp = z(i)
    f_x(i + 1) = (sigma * (y_temp - x_temp)) * dt
    f_y(i + 1) = (x_temp * (rho - z_temp) - y_temp) * dt
    f_z(i + 1) = (x_temp * y_temp - beta * z_temp) * dt
    k1x = f_x(i + 1)
    k1y = f_y(i + 1)
    k1z = f_z(i + 1)
    x_temp = x(i) + k1x / 2.0_dp
    y_temp = y(i) + k1y / 2.0_dp
    z_temp = z(i) + k1z / 2.0_dp
    f_x(i + 1) = (sigma * (y_temp - x_temp)) * dt
    f_y(i + 1) = (x_temp * (rho - z_temp) - y_temp) * dt
    f_z(i + 1) = (x_temp * y_temp - beta * z_temp) * dt
    k2x = f_x(i + 1) !Scratch
    k2y = f_y(i + 1)
    k2z = f_z(i + 1)
    x_temp = x(i) + k2x / 2.0_dp
    y_temp = y(i) + k2y / 2.0_dp
    z_temp = z(i) + k2z / 2.0_dp
    f_x(i + 1) = (sigma * (y_temp - x_temp)) * dt
    f_y(i + 1) = (x_temp * (rho - z_temp) - y_temp) * dt
    f_z(i + 1) = (x_temp * y_temp - beta * z_temp) * dt
    k3x = f_x(i + 1)
    k3y = f_y(i + 1)
    k3z = f_z(i + 1)
    x_temp = x(i) + k3x
    y_temp = y(i) + k3y
    z_temp = z(i) + k3z
    f_x(i + 1) = (sigma * (y_temp - x_temp)) * dt
    f_y(i + 1) = (x_temp * (rho - z_temp) - y_temp) * dt
    f_z(i + 1) = (x_temp * y_temp - beta * z_temp) * dt
    k4x = f_x(i + 1)
    k4y = f_y(i + 1)
    k4z = f_z(i + 1)
    IF (i > 1) THEN
      IF (f_z(i - 1) > 0 .AND. f_z(i) < 0) THEN
        z_m(j) = ABS((z(i - 1) + z(i)) / 2.0_dp)
        IF (j >= 3) WRITE(max_dat, my_format) z_m(j - 1), z_m(j)
        j = j + 1
      END IF
    END IF
    x(i + 1) = x(i) + (k1x + 2.0_dp * k2x + 2.0_dp * k3x + k4x) / 6.0_dp
    y(i + 1) = y(i) + (k1y + 2.0_dp * k2y + 2.0_dp * k3y + k4y) / 6.0_dp
    z(i + 1) = z(i) + (k1z + 2.0_dp * k2z + 2.0_dp * k3z + k4z) / 6.0_dp
    i = i + 1
  END DO lorenz
  CLOSE(max_dat)
  DEALLOCATE(f_x, x, f_y, y, f_z, z, z_m)
END PROGRAM discrete