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, r, r_read, b, b_read, dt, dt_read,    &
    x_read, x_temp, y_read, y_temp, z_read, z_temp
  REAL(dp), DIMENSION(:), ALLOCATABLE, ASYNCHRONOUS :: f_x, x, f_y, y, f_z, z, z_m
  REAL(dp), PARAMETER :: SIX = 6.0
  CHARACTER(*), PARAMETER :: my_format = '(2f12.8)'
END MODULE Shared
!1. dp?, D+00?
!2. i + 1?
!3. Stupid things with _plot?
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))
  WRITE(6, '(a4)', advance='no') "x = "
  READ *, x_read
  x(1) = x_read
  WRITE(6, '(a4)', advance='no') "y = "
  READ *, y_read
  y(1) = y_read
  WRITE(6, '(a4)', advance='no') "z = "
  READ *, z_read
  z(1) = z_read
  WRITE(6, '(a8)', advance='no') "sigma = "
  READ *, sigma_read
  sigma = sigma_read
  WRITE(6, '(a4)', advance='no') "r = "
  READ *, r_read
  r = r_read
  WRITE(6, '(a4)', advance='no') "b = "
  READ *, b_read
  b = b_read
  WRITE(6, '(a5)', advance='no') "dt = "
  READ *, dt_read
  dt = dt_read
  i = 1
  j = 2
  OPEN(max_dat, file='max.dat', access='sequential', status='unknown', action='write', asynchronous='yes')
  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 * (r - z_temp) - y_temp) * dt
    f_z(i + 1) = (x_temp * y_temp - b * 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 * (r - z_temp) - y_temp) * dt
    f_z(i + 1) = (x_temp * y_temp - b * z_temp) * dt
    k2x = f_x(i + 1)
    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 * (r - z_temp) - y_temp) * dt
    f_z(i + 1) = (x_temp * y_temp - b * 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 * (r - z_temp) - y_temp) * dt
    f_z(i + 1) = (x_temp * y_temp - b * 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, asynchronous='yes') 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) / SIX
    y(i + 1) = y(i) + (k1y + 2.0_dp * k2y + 2.0_dp * k3y + k4y) / SIX
    z(i + 1) = z(i) + (k1z + 2.0_dp * k2z + 2.0_dp * k3z + k4z) / SIX
    i = i + 1
  END DO lorenz
  CLOSE(max_dat)
  DEALLOCATE(f_x, x, f_y, y, f_z, z, z_m)
END PROGRAM discrete