program main

!*****************************************************************************80
!
!! lorenz96_ode() solves a version of the Lorenz96 ODE.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    30 September 2025
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  integer, parameter :: n = 4
  integer, parameter :: nt = 2000

  character ( len = 255 ) command_filename
  integer command_unit
  character ( len = 255 ) data_filename
  integer data_unit
  real ( kind = rk8 ) dt
  real ( kind = rk8 ), parameter :: force = 8.0D+00
  integer i
  integer j
  external lorenz96_rhs
  real ( kind = rk8 ), parameter :: perturb = 0.001D+00
  real ( kind = rk8 ) r8_normal_01
  real ( kind = rk8 ) t(0:nt)
  real ( kind = rk8 ), parameter :: t0 = 0.0D+00
  real ( kind = rk8 ), parameter :: tstop = 30.0D+00
  real ( kind = rk8 ) y(n,0:nt)

  call timestamp ( )
  write ( *, '(a)' ) ''
  write ( *, '(a)' ) 'lorenz96_ode():'
  write ( *, '(a)' ) '  Fortran90 version'
  write ( *, '(a)' ) '  Compute solutions of the Lorenz96 system.'
  write ( *, '(a)' ) '  Write data to a file for use by gnuplot().'
!
!  Data
!
  dt = ( tstop - t0 ) / real ( nt, kind = rk8 )
!
!  Initial conditions.
!
  do j = 0, nt
    t(j) = &
      ( real ( nt - j, kind = rk8 ) * t0 &
      + real (      j, kind = rk8 ) * tstop ) &
      / real ( nt,     kind = rk8 )
  end do

  do i = 1, n
    y(i,0) = force + perturb * r8_normal_01 ( )
  end do
!
!  Compute the approximate solution at equally spaced times.
!
  do j = 0, nt - 1

    call rk4vec ( t(j), n, y(1:n,j), dt, lorenz96_rhs, y(1:n,j+1) )

  end do
!
!  Create the data file.
!
  call get_unit ( data_unit )
  data_filename = 'lorenz96_ode_data.txt'
  open ( unit = data_unit, file = data_filename, status = 'replace' )
  do j = 0, nt
    write ( data_unit, '(2x,g14.6,2x,g14.6,2x,g14.6,2x,g14.6,2x,g14.6)' ) t(j), y(1:n,j)
  end do
  close ( unit = data_unit )
  write ( *, '(a)' ) '  Created data file "' // trim ( data_filename ) // '".'
!
!  Create the command file.
!
  call get_unit ( command_unit )
  command_filename = 'lorenz96_ode_commands.txt'
  open ( unit = command_unit, file = command_filename, status = 'replace' )
  write ( command_unit, '(a)' ) '# ' // trim ( command_filename )
  write ( command_unit, '(a)' ) '#'
  write ( command_unit, '(a)' ) '# Usage:'
  write ( command_unit, '(a)' ) '#  gnuplot < ' // trim ( command_filename )
  write ( command_unit, '(a)' ) '#'
  write ( command_unit, '(a)' ) 'set term png'
  write ( command_unit, '(a)' ) &
    'set output "lorenz96_time.png"'
  write ( command_unit, '(a)' ) 'set xlabel "<--- T --->"'
  write ( command_unit, '(a)' ) 'set ylabel "<--- Y1(T), Y2(T), Y3(T) --->"'
  write ( command_unit, '(a)' ) &
    'set title "Y1(T), Y2(T), Y3(T) versus Time"'
  write ( command_unit, '(a)' ) 'set grid'
  write ( command_unit, '(a)' ) 'set style data lines'
  write ( command_unit, '(a)' ) 'plot "' // trim ( data_filename ) // &
    '" using 1:2 lw 3 linecolor rgb "blue",' // &
    ' "" using 1:3 lw 3 linecolor rgb "red",' // &
    ' "" using 1:4 lw 3 linecolor rgb "green"'
  write ( command_unit, '(a)' ) &
    'set output "lorenz96_3d.png"'
  write ( command_unit, '(a)' ) 'set xlabel "<--- Y1(T) --->"'
  write ( command_unit, '(a)' ) 'set ylabel "<--- Y2(T) --->"'
  write ( command_unit, '(a)' ) 'set zlabel "<--- Y3(T) --->"'
  write ( command_unit, '(a)' ) &
    'set title "(Y1(T),Y2(T),Y3(T)) trajectory"'
  write ( command_unit, '(a)' ) 'set grid'
  write ( command_unit, '(a)' ) 'set style data lines'
  write ( command_unit, '(a)' ) 'splot "' // trim ( data_filename ) // &
    '" using 2:3:4 lw 1 linecolor rgb "blue"'
  close ( unit = command_unit )
  write ( *, '(a)' ) &
    '  Created command file "' // trim ( command_filename ) // '".'
!
!  Terminate.
!
  write ( *, '(a)' ) ''
  write ( *, '(a)' ) 'lorenz96_ode():'
  write ( *, '(a)' ) '  Normal end of execution.'
  write ( *, '(a)' ) ''
  call timestamp ( )

  stop
end
subroutine get_unit ( iunit )

!*****************************************************************************80
!
!! get_unit() returns a free Fortran unit number.
!
!  Discussion:
!
!    A "free" Fortran unit number is a value between 1 and 99 which
!    is not currently associated with an I/O device.  A free Fortran unit
!    number is needed in order to open a file with the OPEN command.
!
!    If IUNIT = 0, then no free Fortran unit could be found, although
!    all 99 units were checked (except for units 5, 6 and 9, which
!    are commonly reserved for console I/O).
!
!    Otherwise, IUNIT is a value between 1 and 99, representing a
!    free Fortran unit.  Note that GET_UNIT assumes that units 5 and 6
!    are special, and will never return those values.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    26 October 2008
!
!  Author:
!
!    John Burkardt
!
!  Output:
!
!    integer IUNIT, the free unit number.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  integer i
  integer ios
  integer iunit
  logical lopen

  iunit = 0

  do i = 1, 99

    if ( i /= 5 .and. i /= 6 .and. i /= 9 ) then

      inquire ( unit = i, opened = lopen, iostat = ios )

      if ( ios == 0 ) then
        if ( .not. lopen ) then
          iunit = i
          return
        end if
      end if

    end if

  end do

  return
end
function i4_modp ( i, j )

!*****************************************************************************80
!
!! i4_modp() returns the nonnegative remainder of I4 division.
!
!  Discussion:
!
!    If
!      NREM = I4_MODP ( I, J )
!      NMULT = ( I - NREM ) / J
!    then
!      I = J * NMULT + NREM
!    where NREM is always nonnegative.
!
!    The MOD function computes a result with the same sign as the
!    quantity being divided.  Thus, suppose you had an angle A,
!    and you wanted to ensure that it was between 0 and 360.
!    Then mod(A,360) would do, if A was positive, but if A
!    was negative, your result would be between -360 and 0.
!
!    On the other hand, I4_MODP(A,360) is between 0 and 360, always.
!
!    An I4 is an integer value.
!
!  Example:
!
!        I     J     MOD I4_MODP    Factorization
!
!      107    50       7       7    107 =  2 *  50 + 7
!      107   -50       7       7    107 = -2 * -50 + 7
!     -107    50      -7      43   -107 = -3 *  50 + 43
!     -107   -50      -7      43   -107 =  3 * -50 + 43
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    02 March 1999
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    integer I, the number to be divided.
!
!    integer J, the number that divides I.
!
!  Output:
!
!    integer I4_MODP, the nonnegative remainder when I is divided by J.
!
  implicit none

  integer i
  integer i4_modp
  integer j
  integer value

  if ( j == 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'i4_modp(): Fatal error!'
    write ( *, '(a,i8)' ) '  Illegal divisor J = ', j
    stop 1
  end if

  value = mod ( i, j )

  if ( value < 0 ) then
    value = value + abs ( j )
  end if

  i4_modp = value

  return
end
function i4_wrap ( value, ilo, ihi )

!*****************************************************************************80
!
!! i4_wrap() forces an I4 to lie between given limits by wrapping.
!
!  Discussion:
!
!    An I4 is an integer value.
!
!    There appears to be a bug in the GFortran compiler which can lead to
!    erroneous results when the first argument of i4_wrap() is an expression.
!    In particular:
!
!    do i = 1, 3
!      if ( test ) then
!        i4 = i4_wrap ( i + 1, 1, 3 )
!      end if
!    end do
!
!    was, when I = 3, returning I4 = 3.  So I had to replace this with
!
!    do i = 1, 3
!      if ( test ) then
!        i4 = i + 1
!        i4 = i4_wrap ( i4, 1, 3 )
!      end if
!    end do
!
!  Example:
!
!    ILO = 4, IHI = 8
!
!    I  Value
!
!    -2     8
!    -1     4
!     0     5
!     1     6
!     2     7
!     3     8
!     4     4
!     5     5
!     6     6
!     7     7
!     8     8
!     9     4
!    10     5
!    11     6
!    12     7
!    13     8
!    14     4
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    30 September 2025
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    integer value: a value.
!
!    integer ILO, IHI, the desired bounds.
!
!  Output:
!
!    integer I4_WRAP, a "wrapped" version of the value.
!
  implicit none

  integer i4_modp
  integer i4_wrap
  integer ihi
  integer ilo
  integer value

  i4_wrap = ilo + i4_modp ( value - ilo, ihi + 1 - ilo )

  return
end
subroutine lorenz96_rhs ( t, n, y, dydt )

!*****************************************************************************80
!
!! lorenz96_rhs() evaluates the right hand side of the Lorenz96 ODE.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    30 September 2025
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    real ( kind = rk8 ) T, the value of the independent variable.
!
!    integer N, the spatial dimension.
!
!    real ( kind = rk8 ) Y(N), the values of the dependent variables
!    at time T.
!
!  Output:
!
!    real ( kind = rk8 ) DYDT(N), the values of the derivatives
!    of the dependent variables at time T.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  integer n

  real ( kind = rk8 ) dydt(n)
  real ( kind = rk8 ), parameter :: force = 8.0D+00
  integer i
  integer i4_wrap
  integer im1
  integer im2
  integer ip1
  real ( kind = rk8 ) t
  real ( kind = rk8 ) y(n)

  call r8_fake_use ( t )

  do i = 1, n
    ip1 = i4_wrap ( i + 1, 1, n )
    im1 = i4_wrap ( i - 1, 1, n )
    im2 = i4_wrap ( i - 2, 1, n )
    dydt(i) = ( y(ip1) - y(im2) ) * y(im1) - y(i) + force
  end do

  return
end
subroutine r8_fake_use ( x )

!*****************************************************************************80
!
!! r8_fake_use() pretends to use an R8 variable.
!
!  Discussion:
!
!    Some compilers will issue a warning if a variable is unused.
!    Sometimes there's a good reason to include a variable in a program,
!    but not to use it.  Calling this function with that variable as
!    the argument will shut the compiler up.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    21 April 2020
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    real ( kind = rk8 ) X, the variable to be "used".
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  real ( kind = rk8 ) x

  if ( x /= x ) then
    write ( *, '(a)' ) '  r8_fake_use(): variable is NAN.'
  end if

  return
end
function r8_normal_01 ( )

!*****************************************************************************80
!
!! r8_normal_01() returns a unit pseudonormal R8.
!
!  Discussion:
!
!    The standard normal probability distribution function (PDF) has
!    mean 0 and standard deviation 1.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    06 August 2013
!
!  Author:
!
!    John Burkardt
!
!  Output:
!
!    real ( kind = rk8 ) R8_NORMAL_01, a sample of the standard
!    normal PDF.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  real ( kind = rk8 ) r1
  real ( kind = rk8 ) r2
  real ( kind = rk8 ) r8_normal_01
  real ( kind = rk8 ), parameter :: r8_pi = 3.141592653589793D+00
  real ( kind = rk8 ) x

  call random_number ( harvest = r1 )
  call random_number ( harvest = r2 )

  x = sqrt ( - 2.0D+00 * log ( r1 ) ) * cos ( 2.0D+00 * r8_pi * r2 )

  r8_normal_01 = x

  return
end
subroutine rk4vec ( t0, m, u0, dt, f, u )

!*****************************************************************************80
!
!! rk4vec() takes one Runge-Kutta step for a vector system.
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    08 October 2013
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    real ( kind = rk8 ) T0, the current time.
!
!    integer M, the spatial dimension
!
!    real ( kind = rk8 ) U0(N), the solution estimate at the current time.
!
!    real ( kind = rk8 ) DT, the time step.
!
!    external F, a subroutine of the form 
!      subroutine f ( t, m, u, uprime ) 
!    which evaluates the derivative UPRIME(1:N) given the time T and
!    solution vector U(1:N).
!
!  Output:
!
!    real ( kind = rk8 ) U(M), the fourth-order Runge-Kutta solution 
!    estimate at time T0+DT.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  integer m

  real ( kind = rk8 ) dt
  external f
  real ( kind = rk8 ) f0(m)
  real ( kind = rk8 ) f1(m)
  real ( kind = rk8 ) f2(m)
  real ( kind = rk8 ) f3(m)
  real ( kind = rk8 ) t0
  real ( kind = rk8 ) t1
  real ( kind = rk8 ) t2
  real ( kind = rk8 ) t3
  real ( kind = rk8 ) u(m)
  real ( kind = rk8 ) u0(m)
  real ( kind = rk8 ) u1(m)
  real ( kind = rk8 ) u2(m)
  real ( kind = rk8 ) u3(m)
!
!  Get four sample values of the derivative.
!
  call f ( t0, m, u0, f0 )

  t1 = t0 + dt / 2.0D+00
  u1(1:m) = u0(1:m) + dt * f0(1:m) / 2.0D+00
  call f ( t1, m, u1, f1 )

  t2 = t0 + dt / 2.0D+00
  u2(1:m) = u0(1:m) + dt * f1(1:m) / 2.0D+00
  call f ( t2, m, u2, f2 )

  t3 = t0 + dt
  u3(1:m) = u0(1:m) + dt * f2(1:m)
  call f ( t1, m, u1, f3 )
!
!  Combine them to estimate the solution at time T0 + DT.
!
  u(1:m) = u0(1:m) + dt * ( f0(1:m) + 2.0D+00 * f1(1:m) + 2.0D+00 * f2(1:m) &
    + f3(1:m) ) / 6.0D+00

  return
end
subroutine timestamp ( )

!*****************************************************************************80
!
!! timestamp() prints the current YMDHMS date as a time stamp.
!
!  Example:
!
!    31 May 2001   9:45:54.872 AM
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    18 May 2013
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  character ( len = 8 ) ampm
  integer d
  integer h
  integer m
  integer mm
  character ( len = 9 ), parameter, dimension(12) :: month = (/ &
    'January  ', 'February ', 'March    ', 'April    ', &
    'May      ', 'June     ', 'July     ', 'August   ', &
    'September', 'October  ', 'November ', 'December ' /)
  integer n
  integer s
  integer values(8)
  integer y

  call date_and_time ( values = values )

  y = values(1)
  m = values(2)
  d = values(3)
  h = values(5)
  n = values(6)
  s = values(7)
  mm = values(8)

  if ( h < 12 ) then
    ampm = 'AM'
  else if ( h == 12 ) then
    if ( n == 0 .and. s == 0 ) then
      ampm = 'Noon'
    else
      ampm = 'PM'
    end if
  else
    h = h - 12
    if ( h < 12 ) then
      ampm = 'PM'
    else if ( h == 12 ) then
      if ( n == 0 .and. s == 0 ) then
        ampm = 'Midnight'
      else
        ampm = 'AM'
      end if
    end if
  end if

  write ( *, '(i2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) &
    d, trim ( month(m) ), y, h, ':', n, ':', s, '.', mm, trim ( ampm )

  return
end