Fortran Aide-mémoire

program hello
  ! A complete Fortran program structure
  implicit none
  ! declarations go here
  integer :: status = 0

  print *, "Hello, World!"   ! list-directed output to stdout

  ! executable statements
  call do_work(status)

  print *, "Exit status: ", status
contains
  subroutine do_work(st)
    integer, intent(out) :: st
    st = 0
    print *, "Working..."
  end subroutine do_work
end program hello

! Compile: gfortran hello.f90 -o hello
! Run:     ./hello

program variables
  implicit none
  ! Integer types
  integer :: count = 0
  integer(kind=8) :: big = 9223372036854775807_8   ! 64-bit
  ! Real types
  real :: x = 3.14              ! default (often 32-bit)
  real(kind=8) :: y = 2.718281828459045_8   ! double precision
  double precision :: z = 1.0d0
  ! Complex
  complex :: c = (1.0, 2.0)    ! 1 + 2i
  complex(kind=8) :: cw = (1.0_8, 2.0_8)
  ! Character
  character(len=20) :: name = "Alice"
  character(len=:), allocatable :: flexible   ! deferred length
  ! Logical
  logical :: flag = .true.
  ! Print all
  print *, count, big
  print *, x, y, z
  print *, c, cw
  print *, name, flag
end program variables

program constants
  implicit none
  ! Named constants via 'parameter' attribute
  integer, parameter :: MAX_SIZE = 100
  real, parameter :: PI = 3.14159265
  real, parameter :: E = 2.718281828
  character(*), parameter :: APP_NAME = "MyApp"
  ! Using parameters
  real :: arr(MAX_SIZE)
  arr = 0.0
  print *, APP_NAME, " size=", MAX_SIZE
  print *, "Circumference: ", 2.0 * PI * 5.0

  ! Enum-like via parameter
  integer, parameter :: SUNDAY = 1, MONDAY = 2, TUESDAY = 3
  integer :: day = MONDAY
  print *, "Day code: ", day
end program constants

program operators
  implicit none
  integer :: a = 17, b = 5
  real :: x = 2.0
  ! Arithmetic
  print *, a + b, a - b, a * b    ! 22 12 85
  print *, a / b                   ! 3 (integer division!)
  print *, real(a) / b             ! 3.4 (cast to real)
  print *, a ** 2                  ! 289 (exponentiation)
  print *, mod(a, b)               ! 2 (modulo)
  print *, modulo(a, b)            ! 2 (differs for negatives)
  ! Relational (both forms work)
  print *, a > b, a < b            ! T F
  print *, a .gt. b, a .lt. b      ! T F (old form)
  print *, a == b, a /= b          ! F T
  ! Logical
  print *, (a > 0) .and. (b > 0)   ! T
  print *, (a > 0) .or. (b < 0)    ! T
  print *, .not. (a > 0)           ! F
  print *, (a > 0) .eqv. (b > 0)   ! T (equivalence)
  ! String concatenation
  character(10) :: s1 = "Hello", s2 = "World"
  print *, trim(s1) // " " // trim(s2)   ! Hello World
end program operators

program intrinsics
  implicit none
  real :: x = -3.7, y = 2.5
  ! Math functions
  print *, abs(x)         ! 3.7
  print *, sqrt(2.0)      ! 1.414...
  print *, exp(1.0)       ! 2.718... (e^x)
  print *, log(2.0)       ! 0.693... (natural log)
  print *, log10(1000.0)  ! 3.0
  print *, sin(3.14159/2) ! 1.0
  print *, cos(0.0)       ! 1.0
  print *, atan2(1.0,1.0) ! 0.785... (pi/4)
  ! Rounding
  print *, int(x)         ! -3 (truncate toward zero)
  print *, nint(x)        ! -4 (nearest integer)
  print *, floor(x)       ! -4 (toward -inf)
  print *, ceiling(x)     ! -3 (toward +inf)
  print *, abs(x), max(x, y), min(x, y)   ! 3.7 2.5 -3.7
  ! Type conversion
  print *, real(5)        ! 5.0
  print *, int(3.9)       ! 3
  ! Inquiry
  real :: arr(10)
  print *, size(arr)      ! 10
  print *, huge(1)        ! 2147483647
  print *, tiny(1.0)      ! smallest positive real
end program intrinsics

program if_demo
  implicit none
  integer :: score = 85
  character(1) :: grade

  ! Multi-branch if/else if/else
  if (score >= 90) then
    grade = 'A'
  else if (score >= 80) then
    grade = 'B'
  else if (score >= 70) then
    grade = 'C'
  else if (score >= 60) then
    grade = 'D'
  else
    grade = 'F'
  end if

  print *, "Score ", score, " -> Grade ", grade

  ! Logical if (single statement, no 'then')
  if (score < 0 .or. score > 100) stop "Invalid score"

  ! Arithmetic if (OBSOLETE - avoid)
  ! if (x) 10, 20, 30   ! jump to label based on sign
end program if_demo

program case_demo
  implicit none
  integer :: day = 3
  character(1) :: ch = 'A'
  character(10) :: day_name

  ! Integer select case
  select case (day)
  case (1)
    day_name = "Monday"
  case (2)
    day_name = "Tuesday"
  case (3)
    day_name = "Wednesday"
  case (4:5)
    day_name = "Thu/Fri"
  case (6:7)
    day_name = "Weekend"
  case default
    day_name = "Invalid"
  end select
  print *, day_name

  ! Character select case (case-insensitive via pre-upper)
  select case (ch)
  case ('A':'Z')
    print *, "Uppercase letter"
  case ('a':'z')
    print *, "Lowercase letter"
  case ('0':'9')
    print *, "Digit"
  case default
    print *, "Other"
  end select

  ! Logical select case
  select case (day > 5)
  case (.true.)
    print *, "Weekend!"
  case (.false.)
    print *, "Weekday"
  end select
end program case_demo

program do_loops
  implicit none
  integer :: i, j, total

  ! Basic counted loop: do var = start, end [, step]
  do i = 1, 5
    print *, i          ! 1 2 3 4 5
  end do

  ! With step
  do i = 10, 1, -1      ! countdown
    print *, i
  end do

  do i = 0, 100, 25     ! 0 25 50 75 100
    print *, i
  end do

  ! Implied-do (inline, for array init / I/O)
  integer :: arr(5) = [(i**2, i=1,5)]   ! 1 4 9 16 25
  print *, arr
  print *, (i, i=1,3)                   ! 1 2 3

  ! Nested loops with labels (for cycle/exit targeting)
  total = 0
  outer: do i = 1, 3
    inner: do j = 1, 3
      total = total + i*j
    end do inner
  end do outer
  print *, "Total: ", total
end program do_loops

program while_demo
  implicit none
  integer :: n, count
  real :: x, sum

  ! do while: pre-test loop
  n = 1024
  count = 0
  do while (n > 1)
    n = n / 2
    count = count + 1
  end do
  print *, "log2(1024) = ", count   ! 10

  ! Infinite loop with exit
  sum = 0.0
  do
    read(*, *) x
    if (x < 0) exit          ! leave loop
    sum = sum + x
  end do
  print *, "Sum: ", sum

  ! do ... end do with conditional exit
  n = 1
  factorial_loop: do
    if (n > 10) exit factorial_loop
    print *, n, factorial(n)
    n = n + 1
  end do factorial_loop

contains
  recursive function factorial(n) result(f)
    integer, intent(in) :: n
    integer :: f
    if (n <= 1) then
      f = 1
    else
      f = n * factorial(n-1)
    end if
  end function factorial
end program while_demo

program loop_control
  implicit none
  integer :: i, j

  ! cycle: skip to next iteration (like 'continue' in C)
  do i = 1, 10
    if (mod(i, 2) == 0) cycle   ! skip even numbers
    print *, i                  ! 1 3 5 7 9
  end do

  ! exit: break out of loop (like 'break' in C)
  do i = 1, 100
    if (i * i > 50) then
      print *, "Stopped at i=", i
      exit
    end if
  end do

  ! Named loops: cycle/exit can target outer loops
  search: do i = 1, 5
    do j = 1, 5
      if (i + j == 7) then
        print *, "Found: ", i, "+", j, "= 7"
        exit search             ! break out of OUTER loop
      end if
    end do
  end do search

  ! Early exit from a search
  integer :: arr(10) = [3, 1, 4, 1, 5, 9, 2, 6, 5, 3]
  do i = 1, size(arr)
    if (arr(i) == 9) then
      print *, "Found 9 at index ", i
      exit
    end if
  end do
end program loop_control

program array_decl
  implicit none
  ! Declaration with dimension
  integer :: a(5)                  ! 1D, indices 1..5
  integer :: b(0:4)                ! 1D, indices 0..4 (custom lower bound)
  real :: c(3, 4)                  ! 2D, 3 rows x 4 cols
  real, dimension(10) :: d         ! using dimension attribute

  ! Initialization at declaration
  integer :: x(5) = [1, 2, 3, 4, 5]
  integer :: y(5) = [(i*2, i=1,5)] ! implied-do: 2 4 6 8 10
  integer :: z(5) = 0              ! all zeros
  real :: m(2,2) = reshape([1,2,3,4], [2,2])

  ! Allocation later
  print *, size(x), lbound(x), ubound(x)   ! 5 1 5
  print *, size(c, dim=1)                  ! 3 (rows)
  print *, size(c, dim=2)                  ! 4 (cols)

  ! Array of characters
  character(10) :: names(3) = ["Alice", "Bob", "Carol"]
  print *, names(2)                        ! Bob
end program array_decl

program array_sections
  implicit none
  integer :: a(10) = [(i, i=1,10)]
  integer :: b(5)
  integer :: idx(3) = [2, 5, 7]

  ! Array sections (slicing): a(start:end[:stride])
  print *, a(3:7)         ! 3 4 5 6 7
  print *, a(1:10:2)      ! 1 3 5 7 9 (stride 2)
  print *, a(:5)          ! 1 2 3 4 5 (start defaults to 1)
  print *, a(6:)          ! 6 7 8 9 10 (end defaults to ubound)
  print *, a(::3)         ! 1 4 7 10
  print *, a(10:1:-1)     ! 10 9 8 ... 1 (reverse)

  ! Vector subscript (gather)
  b = a(idx)              ! b = [a(2), a(5), a(7)] = [2, 5, 7]
  print *, b

  ! Assigning to a section
  a(2:4) = [99, 98, 97]
  print *, a(1:5)         ! 1 99 98 97 5

  ! 2D sections
  real :: m(4,4) = 0.0
  m(2:3, 2:3) = 1.0       ! set 2x2 sub-block
  print *, m(2,2), m(3,3) ! 1.0 1.0

  ! Where construct (masked array assignment)
  where (a > 50) a = 0    ! set elements > 50 to 0
end program array_sections

program matrices
  implicit none
  real :: A(3,3), B(3,3), C(3,3)
  integer :: i, j

  ! Initialize with implied-do
  A = reshape([(real(i), i=1,9)], [3,3])
  ! A = 1 4 7
  !     2 5 8
  !     3 6 9   (column-major fill!)

  ! Identity matrix
  B = 0.0
  do i = 1, 3
    B(i,i) = 1.0
  end do

  ! Matrix multiplication (intrinsic)
  C = matmul(A, B)        ! A * I = A
  print *, C(1,1), C(2,2) ! 1.0 5.0

  ! Transpose
  print *, transpose(A)(1,:)   ! 1 2 3 (first row of A^T)

  ! Element-wise operations
  C = A + B               ! element-wise add
  C = A * 2.0             ! scalar multiply
  C = A * B               ! element-wise (NOT matmul!)

  ! Array reduction along a dimension
  print *, sum(A, dim=1)   ! column sums: 6 15 24
  print *, sum(A, dim=2)   ! row sums: 12 15 18
  print *, maxval(A)       ! 9.0
  print *, maxloc(A)       ! 3 3 (location of max)

  ! Reshape
  integer :: flat(6) = [1,2,3,4,5,6]
  integer :: mat(2,3)
  mat = reshape(flat, [2,3])
end program matrices

program alloc_demo
  implicit none
  integer, allocatable :: arr(:), matrix(:,:)
  integer :: n, m, i, stat

  ! Get size from user
  print *, "Enter size:"
  read(*, *) n
  m = n * 2

  ! Allocate
  allocate(arr(n), matrix(n, m), stat=stat)
  if (stat /= 0) then
    print *, "Allocation failed!"
    stop 1
  end if

  ! Use the arrays
  arr = [(i, i=1, n)]
  matrix = 0.0
  do i = 1, n
    matrix(i, :) = i
  end do

  print *, size(arr), size(matrix, dim=2)
  print *, allocated(arr)   ! T

  ! Deallocate (or let it auto-deallocate at scope exit)
  deallocate(arr, matrix)
  print *, allocated(arr)   ! F

  ! Automatic reallocation on assignment (Fortran 2003)
  integer, allocatable :: flex(:)
  flex = [1, 2, 3]          ! auto-allocates to size 3
  flex = [flex, 4, 5]       ! reallocates to size 5: 1 2 3 4 5
  print *, flex
  deallocate(flex)
end program alloc_demo

program array_funcs
  implicit none
  integer :: a(5) = [3, 1, 4, 1, 5, 9, 2, 6]
  ! Wait, let's fix size
  integer :: b(8) = [3, 1, 4, 1, 5, 9, 2, 6]
  integer :: c(5) = [10, 20, 30, 40, 50]
  logical :: mask(5) = [.true., .false., .true., .false., .true.]

  ! Inquiry
  print *, size(b)        ! 8
  print *, shape(b)       ! 8
  print *, lbound(b), ubound(b)   ! 1 8

  ! Reductions
  print *, sum(b)         ! 31
  print *, product(c)     ! 120000000
  print *, maxval(b), minval(b)   ! 9 1
  print *, maxloc(b)      ! 6 (index of max)
  print *, minloc(b)      ! 2 (index of min, first occurrence)
  print *, count(b > 3)   ! 5 (number of true elements)
  print *, any(b > 8)     ! T (at least one)
  print *, all(b > 0)     ! T (all of them)

  ! With mask
  print *, sum(b, mask=b > 3)     ! sum of elements > 3
  print *, pack(b, b > 3)         ! compact array of elements > 3
  print *, unpack([1,2], mask, 0) ! spread values per mask

  ! Manipulation
  print *, cshift(b, 2)   ! circular shift left by 2
  print *, eoshift(b, 2)  ! end-off shift left by 2 (fills 0)
  print *, merge(b, c, mask)  ! element-wise: b where mask true, else c

  ! Sorting (Fortran 2003+)
  integer :: sorted(8)
  sorted = b
  call sort_array(sorted)   ! custom sort (no built-in sort)
end program array_funcs

program char_decl
  implicit none
  ! Fixed-length strings
  character(10) :: s1 = "Hello"
  character(len=20) :: s2 = "World"
  character(20) :: s3          ! len= keyword optional

  ! Deferred-length (allocatable) - Fortran 2003+
  character(:), allocatable :: flex
  flex = "Dynamic"             ! len=7
  flex = "Now longer string"   ! reallocates to len=18
  print *, len(flex)           ! 18

  ! Array of strings
  character(15) :: names(3) = ["Alice", "Bob", "Carol"]

  ! Single character
  character(1) :: ch = 'A'
  character :: ch2 = 'B'       ! len=1 default

  print *, s1                  ! "Hello     " (padded to 10)
  print *, trim(s1)            ! "Hello" (no padding)
  print *, len(s1), len_trim(s1)   ! 10 5
  print *, names(2)            ! "Bob" (padded)
end program char_decl

program string_ops
  implicit none
  character(20) :: first = "John", last = "Doe"
  character(50) :: full

  ! Concatenation with //
  full = first // " " // last   ! "John Doe"
  print *, trim(full)

  ! Repeat
  print *, repeat("-", 30)      ! 30 dashes

  ! Substring (1-indexed, inclusive)
  character(20) :: s = "Hello, World!"
  print *, s(1:5)               ! "Hello"
  print *, s(8:12)              ! "World"
  print *, s(8:)                ! "World!"  (to end)
  print *, s(:5)                ! "Hello"   (from start)

  ! Index (find substring)
  print *, index(s, "World")    ! 8 (position, 0 if not found)
  print *, index(s, "world")    ! 0 (case-sensitive)
  print *, scan(s, "aeiou")     ! 2 (first vowel position)
  print *, verify(s, "abcdefg") ! 1 (first char NOT in set)

  ! Length
  print *, len_trim(s)          ! 13
  print *, adjustl(s)           ! left-justify
  print *, adjustr(s)           ! right-justify
end program string_ops

program string_funcs
  implicit none
  character(20) :: s = "Hello World"
  integer :: i
  character(1) :: ch

  ! Character code conversion
  print *, iachar('A')          ! 65 (ASCII code)
  print *, achar(66)            ! 'B' (from ASCII code)
  print *, ichar('A')           ! processor-dependent (use iachar for ASCII)

  ! Case conversion (manual - no built-in)
  do i = 1, len_trim(s)
    ch = s(i:i)
    if (ch >= 'A' .and. ch <= 'Z') then
      s(i:i) = achar(iachar(ch) + 32)   ! to lowercase
    end if
  end do
  print *, s                    ! "hello world"

  ! Comparison
  print *, lge("apple", "banana")   ! F (lexicographic >=)
  print *, lgt("zebra", "apple")    ! T
  print *, lle("abc", "abcd")       ! T (<=)
  print *, llt("abc", "abd")        ! T (<)

  ! String to number conversion (internal read)
  character(20) :: num_str = "3.14159"
  real :: pi_val
  read(num_str, *) pi_val
  print *, pi_val * 2           ! 6.28318

  ! Number to string (internal write)
  character(20) :: out_str
  integer :: n = 42
  write(out_str, '(I0)') n      ! I0 = minimal-width integer
  print *, "Number: " // trim(out_str)
end program string_funcs

program format_demo
  implicit none
  integer :: n = 42
  real :: pi = 3.14159265
  real :: big = 1234567.89
  character(20) :: name = "Alice"

  ! Format specifiers
  ! Iw     - integer, width w
  ! Iw.m   - integer, width w, at least m digits
  ! Fw.d   - fixed-point real, width w, d decimals
  ! Ew.d   - exponential, width w, d decimals
  ! Aw     - character, width w
  ! A      - character, default width
  ! nX     - n spaces
  ! /      - newline

  write(*, '(I5)') n            ! "   42"
  write(*, '(I5.3)') n          ! "  042"
  write(*, '(I0)') n            ! "42" (minimal)
  write(*, '(F10.4)') pi        ! "    3.1416"
  write(*, '(E12.4)') big       ! " 0.1235E+07"
  write(*, '(ES12.4)') big      ! " 1.2346E+06" (scientific)
  write(*, '(A10)') name        ! "     Alice"
  write(*, '(A, I3)') "n=", n   ! "n= 42"

  ! Multiple items
  write(*, '(A, I3, A, F8.4)') "n=", n, " pi=", pi

  ! Repeated format: 3I4 = three integers width 4
  write(*, '(3I4)') 1, 2, 3     ! "   1   2   3"

  ! Newline and spacing
  write(*, '("Name: ", A, /, "Age: ", I3)') name, n

  ! List-directed (default format)
  print *, name, n, pi
end program format_demo

program parse_demo
  implicit none
  character(100) :: line = "name=Bob,age=30,city=NYC"
  character(50) :: token
  integer :: pos, start, end

  ! Split by comma
  start = 1
  do
    end = index(line(start:), ",")
    if (end == 0) then
      ! last token
      token = line(start:)
      call process(trim(token))
      exit
    end if
    token = line(start:start+end-2)
    call process(trim(token))
    start = start + end
  end do

contains
  subroutine process(t)
    character(*), intent(in) :: t
    integer :: eq_pos
    eq_pos = index(t, "=")
    if (eq_pos > 0) then
      print *, "Key: ", trim(t(:eq_pos-1)), &
               " Value: ", trim(t(eq_pos+1:))
    end if
  end subroutine process
end program parse_demo

! Output:
! Key: name Value: Bob
! Key: age Value: 30
! Key: city Value: NYC

program func_demo
  implicit none
  ! Function declared in interface or as external
  print *, add(3, 4)         ! 7
  print *, square(5)         ! 25
  print *, distance(0.0, 0.0, 3.0, 4.0)   ! 5.0

contains
  ! Basic function with result clause
  function add(a, b) result(c)
    integer, intent(in) :: a, b
    integer :: c
    c = a + b
  end function add

  ! Function with result named same as function (legacy)
  function square(x) result(y)
    integer, intent(in) :: x
    integer :: y
    y = x * x
  end function square

  ! Real function
  function distance(x1, y1, x2, y2) result(d)
    real, intent(in) :: x1, y1, x2, y2
    real :: d
    d = sqrt((x2-x1)**2 + (y2-y1)**2)
  end function distance

  ! Array-valued function
  function reverse(arr) result(rev)
    integer, intent(in) :: arr(:)
    integer :: rev(size(arr))
    integer :: i, n
    n = size(arr)
    do i = 1, n
      rev(i) = arr(n-i+1)
    end do
  end function reverse
end program func_demo

program sub_demo
  implicit none
  integer :: x = 10, y = 20
  real :: arr(5) = [1.0, 2.0, 3.0, 4.0, 5.0]

  ! Subroutines are called with 'call'
  call swap(x, y)
  print *, x, y              ! 20 10

  call scale_array(arr, 2.0)
  print *, arr               ! 2 4 6 8 10

  call fill_zero(arr)
  print *, arr               ! 0 0 0 0 0

contains
  subroutine swap(a, b)
    integer, intent(inout) :: a, b   ! read AND write
    integer :: tmp
    tmp = a
    a = b
    b = tmp
  end subroutine swap

  subroutine scale_array(a, factor)
    real, intent(inout) :: a(:)
    real, intent(in) :: factor
    a = a * factor           ! whole-array operation
  end subroutine scale_array

  subroutine fill_zero(a)
    real, intent(out) :: a(:)        ! write-only (output)
    a = 0.0
  end subroutine fill_zero
end program sub_demo

program pure_demo
  implicit none
  real :: a(5) = [1.0, 2.0, 3.0, 4.0, 5.0]
  real :: b(5)
  integer :: i

  ! Pure function: no side effects, no I/O
  b = square_arr(a)
  print *, b                 ! 1 4 9 16 25

  ! Elemental function: works on scalars AND arrays automatically
  b = cube(a)                ! applies cube() element-wise
  print *, b                 ! 1 8 27 64 125

  print *, cube(2.0)         ! also works on scalar: 8.0

contains
  ! Pure: no side effects, no I/O, no stop, only pure calls
  pure function square_arr(x) result(y)
    real, intent(in) :: x(:)
    real :: y(size(x))
    y = x * x
  end function square_arr

  ! Elemental: scalar signature, but works on arrays too
  elemental function cube(x) result(y)
    real, intent(in) :: x
    real :: y
    y = x * x * x
  end function cube

  ! Pure elemental: both
  pure elemental double precision function sq(x) result(y)
    double precision, intent(in) :: x
    double precision :: y
    y = x * x
  end function sq
end program pure_demo

program optional_demo
  implicit none
  ! All arguments after the first optional must also be optional
  print *, greet("Alice")                    ! Hello, Alice!
  print *, greet("Bob", "Hi")                ! Hi, Bob!
  print *, greet("Carol", greeting="Hey")    ! keyword argument
  print *, greet(greeting="Welcome", name="Dave")  ! all keywords

  ! With present() check
  call log_msg("Starting up")
  call log_msg("Error!", level=2)
  call log_msg("Debug info", level=0, file="debug.log")

contains
  function greet(name, greeting) result(msg)
    character(*), intent(in) :: name
    character(*), intent(in), optional :: greeting
    character(50) :: msg
    character(20) :: g

    if (present(greeting)) then
      g = greeting
    else
      g = "Hello"
    end if
    msg = trim(g) // ", " // name // "!"
  end function greet

  subroutine log_msg(message, level, file)
    character(*), intent(in) :: message
    integer, intent(in), optional :: level
    character(*), intent(in), optional :: file
    integer :: lvl
    lvl = 1
    if (present(level)) lvl = level
    print *, "[L", lvl, "] ", trim(message)
  end subroutine log_msg
end program optional_demo

program nested_demo
  implicit none
  integer :: counter = 0   ! host variable

  ! Internal procedures (in contains) access host variables
  call increment
  call increment
  print *, counter          ! 2

  ! Recursive function
  print *, factorial(5)     ! 120
  print *, fib(10)          ! 55

  ! Mutual recursion (needs forward declaration)
  print *, is_even(4)       ! T

contains
  subroutine increment
    counter = counter + 1   ! modifies host's counter
  end subroutine increment

  recursive function factorial(n) result(f)
    integer, intent(in) :: n
    integer :: f
    if (n <= 1) then
      f = 1
    else
      f = n * factorial(n-1)
    end if
  end function factorial

  recursive function fib(n) result(f)
    integer, intent(in) :: n
    integer :: f
    if (n < 2) then
      f = n
    else
      f = fib(n-1) + fib(n-2)
    end if
  end function fib

  ! Mutual recursion with interface
  recursive function is_even(n) result(r)
    integer, intent(in) :: n
    logical :: r
    interface
      recursive function is_odd(m) result(ro)
        integer, intent(in) :: m
        logical :: ro
      end function is_odd
    end interface
    if (n == 0) then
      r = .true.
    else
      r = is_odd(n-1)
    end if
  end function is_even
end program nested_demo

! geometry.f90 - module file
module geometry
  implicit none
  private                      ! default: everything private
  public :: circle_area, circle_perimeter, PI

  ! Module-level constants (persistent)
  real, parameter :: PI = 3.14159265

contains
  function circle_area(r) result(a)
    real, intent(in) :: r
    real :: a
    a = PI * r * r
  end function circle_area

  function circle_perimeter(r) result(p)
    real, intent(in) :: r
    real :: p
    p = 2.0 * PI * r
  end function circle_perimeter
end module geometry

! main.f90 - using the module
program use_module
  use geometry, only: circle_area, circle_perimeter, PI
  implicit none
  real :: r = 5.0
  print *, "Area: ", circle_area(r)         ! 78.5398
  print *, "Perimeter: ", circle_perimeter(r) ! 31.4159
  print *, "PI: ", PI
end program use_module

! Compile: gfortran geometry.f90 main.f90 -o main

module bank_account
  implicit none
  private                       ! default: everything private
  ! Explicitly export:
  public :: account_t, deposit, withdraw, get_balance

  ! Derived type - can expose type but hide internals
  type :: account_t
    private                    ! components are private
    real :: balance = 0.0
    integer :: id = 0
  contains
    procedure :: balance => get_bal   ! type-bound procedure
    final :: cleanup                  ! destructor
  end type account_t

  ! Module-level counter (private, not exported)
  integer, save :: next_id = 1000

contains
  ! Constructor (factory function)
  function create_account(initial) result(acc)
    type(account_t) :: acc
    real, intent(in) :: initial
    acc%balance = initial
    acc%id = next_id
    next_id = next_id + 1
  end function create_account

  subroutine deposit(acc, amount)
    type(account_t), intent(inout) :: acc
    real, intent(in) :: amount
    acc%balance = acc%balance + amount
  end subroutine deposit

  function get_bal(acc) result(b)
    class(account_t), intent(in) :: acc
    real :: b
    b = acc%balance
  end function get_bal

  subroutine cleanup(acc)
    type(account_t) :: acc
    ! cleanup code (e.g., log closure)
  end subroutine cleanup
end module bank_account

module vector_mod
  implicit none
  private
  public :: vector_t, vector_add, vector_scale

  ! Derived type with type parameters (Fortran 2003)
  type :: vector_t
    real, allocatable :: data(:)
    integer :: length = 0
  contains
    procedure :: norm => vector_norm
    procedure :: print => vector_print
  end type vector_t

  ! Constructor interface (overloaded)
  interface vector_t
    procedure vector_from_array
    procedure vector_from_size
  end interface vector_t

contains
  function vector_from_array(arr) result(v)
    real, intent(in) :: arr(:)
    type(vector_t) :: v
    v%data = arr
    v%length = size(arr)
  end function

  function vector_from_size(n, fill) result(v)
    integer, intent(in) :: n
    real, intent(in) :: fill
    type(vector_t) :: v
    allocate(v%data(n))
    v%data = fill
    v%length = n
  end function

  function vector_norm(self) result(n)
    class(vector_t), intent(in) :: self
    real :: n
    n = sqrt(sum(self%data**2))
  end function

  subroutine vector_print(self)
    class(vector_t), intent(in) :: self
    print *, "Vector(len=", self%length, "): ", self%data
  end subroutine

  function vector_add(a, b) result(c)
    type(vector_t), intent(in) :: a, b
    type(vector_t) :: c
    c%data = a%data + b%data
    c%length = a%length
  end function

  function vector_scale(a, s) result(c)
    type(vector_t), intent(in) :: a
    real, intent(in) :: s
    type(vector_t) :: c
    c%data = a%data * s
    c%length = a%length
  end function
end module vector_mod

module generics_mod
  implicit none
  private
  public :: print_value, add

  ! Generic interface: one name, multiple specific procedures
  interface print_value
    module procedure print_int
    module procedure print_real
    module procedure print_str
    module procedure print_int_array
  end interface print_value

  interface add
    module procedure add_int
    module procedure add_real
    module procedure add_arrays
  end interface add

contains
  subroutine print_int(x)
    integer, intent(in) :: x
    print *, "Integer: ", x
  end subroutine

  subroutine print_real(x)
    real, intent(in) :: x
    print *, "Real: ", x
  end subroutine

  subroutine print_str(s)
    character(*), intent(in) :: s
    print *, "String: ", trim(s)
  end subroutine

  subroutine print_int_array(arr)
    integer, intent(in) :: arr(:)
    print *, "Array: ", arr
  end subroutine

  function add_int(a, b) result(c)
    integer, intent(in) :: a, b
    integer :: c
    c = a + b
  end function

  function add_real(a, b) result(c)
    real, intent(in) :: a, b
    real :: c
    c = a + b
  end function

  function add_arrays(a, b) result(c)
    real, intent(in) :: a(:), b(:)
    real :: c(size(a))
    c = a + b
  end function
end module generics_mod

program use_generics
  use generics_mod
  implicit none
  call print_value(42)              ! Integer: 42
  call print_value(3.14)            ! Real: 3.14
  call print_value("Hello")         ! String: Hello
  call print_value([1,2,3])         ! Array: 1 2 3
  print *, add(2, 3)                ! 5
  print *, add(2.5, 3.5)            ! 6.0
end program use_generics

module vec_ops
  implicit none
  private
  public :: vec3, operator(+), operator(*), assignment(=)

  type :: vec3
    real :: x, y, z
  end type vec3

  ! Overload the + operator for vec3 + vec3
  interface operator(+)
    module procedure vec_add
  end interface

  ! Overload * for vec3 * scalar and scalar * vec3
  interface operator(*)
    module procedure vec_scale_r
    module procedure vec_scale_l
    module procedure vec_dot
  end interface

  ! Overload = for array-to-vec assignment
  interface assignment(=)
    module procedure arr_to_vec
  end interface

contains
  function vec_add(a, b) result(c)
    type(vec3), intent(in) :: a, b
    type(vec3) :: c
    c = vec3(a%x+b%x, a%y+b%y, a%z+b%z)
  end function

  function vec_scale_r(v, s) result(r)
    type(vec3), intent(in) :: v
    real, intent(in) :: s
    type(vec3) :: r
    r = vec3(v%x*s, v%y*s, v%z*s)
  end function

  function vec_scale_l(s, v) result(r)
    real, intent(in) :: s
    type(vec3), intent(in) :: v
    type(vec3) :: r
    r = vec3(v%x*s, v%y*s, v%z*s)
  end function

  ! vec3 * vec3 = dot product (scalar)
  function vec_dot(a, b) result(d)
    type(vec3), intent(in) :: a, b
    real :: d
    d = a%x*b%x + a%y*b%y + a%z*b%z
  end function

  subroutine arr_to_vec(v, arr)
    type(vec3), intent(out) :: v
    real, intent(in) :: arr(3)
    v = vec3(arr(1), arr(2), arr(3))
  end subroutine
end module vec_ops

program use_ops
  use vec_ops
  implicit none
  type(vec3) :: a, b, c
  a = vec3(1.0, 2.0, 3.0)
  b = vec3(4.0, 5.0, 6.0)
  c = a + b              ! vec3 addition
  print *, c%x, c%y, c%z ! 5 7 9
  print *, a * 2.0       ! scale: 2 4 6
  print *, 3.0 * b       ! scale: 12 15 18
  print *, a * b         ! dot product: 32
end program use_ops

program derived_types
  implicit none
  ! Basic derived type (struct)
  type :: Point
    real :: x, y
  end type Point

  ! Type with default initialization
  type :: Person
    character(20) :: name = "Unknown"
    integer :: age = 0
    logical :: active = .true.
  end type Person

  ! Declare and construct
  type(Point) :: p1, p2
  type(Person) :: alice, bob

  ! Structure constructor
  p1 = Point(3.0, 4.0)
  p2 = Point(0.0, 0.0)
  alice = Person("Alice", 30, .true.)
  bob = Person("Bob", 25)        ! uses default for 'active'

  ! Access components with %
  print *, p1%x, p1%y            ! 3.0 4.0
  print *, alice%name, alice%age ! Alice 30
  print *, bob%active            ! T (default)

  ! Modify components
  p1%x = p1%x + 1.0
  alice%age = 31

  ! Whole-type assignment (component-wise copy)
  p2 = p1
  print *, p2%x                  ! 4.0

  ! Array of derived types
  type(Point) :: points(3)
  points(1) = Point(1.0, 1.0)
  points(2) = Point(2.0, 2.0)
  points(3) = Point(3.0, 3.0)
  print *, points(2)%y           ! 2.0
end program derived_types

program type_components
  implicit none
  ! Type with various component kinds
  type :: Student
    integer :: id
    character(20) :: name
    real, allocatable :: grades(:)    ! allocatable component
    integer :: num_grades = 0
  contains
    procedure :: add_grade
    procedure :: average => student_avg
    final :: student_finalize
  end type Student

  ! Overloaded constructor
  interface Student
    module procedure student_init
  end interface Student

  type(Student) :: s

  ! Use custom constructor
  s = Student(101, "Alice")
  call s%add_grade(85.0)
  call s%add_grade(92.0)
  call s%add_grade(78.0)
  print *, s%average()    ! 85.0

contains
  function student_init(id, name) result(s)
    integer, intent(in) :: id
    character(*), intent(in) :: name
    type(Student) :: s
    s%id = id
    s%name = name
    s%num_grades = 0
  end function

  subroutine add_grade(self, g)
    class(Student), intent(inout) :: self
    real, intent(in) :: g
    integer :: n
    n = self%num_grades
    if (n == 0) then
      allocate(self%grades(1))
    else
      ! grow array
      block
        real, allocatable :: tmp(:)
        tmp = self%grades
        deallocate(self%grades)
        allocate(self%grades(n+1))
        self%grades(1:n) = tmp
      end block
    end if
    self%num_grades = n + 1
    self%grades(n+1) = g
  end subroutine

  function student_avg(self) result(avg)
    class(Student), intent(in) :: self
    real :: avg
    if (self%num_grades > 0) then
      avg = sum(self%grades) / self%num_grades
    else
      avg = 0.0
    end if
  end function

  subroutine student_finalize(self)
    type(Student) :: self
    if (allocated(self%grades)) deallocate(self%grades)
  end subroutine
end program type_components

module stack_mod
  implicit none
  private
  public :: stack_t

  type :: stack_t
    integer, allocatable :: data(:)
    integer :: top = 0
    integer :: capacity = 0
  contains
    procedure :: push => stack_push
    procedure :: pop => stack_pop
    procedure :: peek => stack_peek
    procedure :: is_empty => stack_is_empty
    procedure :: size => stack_size
    procedure :: clear => stack_clear
    final :: stack_finalize
  end type stack_t

  interface stack_t
    module procedure stack_init
  end interface

contains
  function stack_init(initial_cap) result(s)
    integer, intent(in), optional :: initial_cap
    type(stack_t) :: s
    integer :: cap
    cap = 16
    if (present(initial_cap)) cap = initial_cap
    allocate(s%data(cap))
    s%capacity = cap
    s%top = 0
  end function

  subroutine stack_push(self, val)
    class(stack_t), intent(inout) :: self
    integer, intent(in) :: val
    if (self%top >= self%capacity) then
      ! grow
      block
        integer, allocatable :: tmp(:)
        tmp = self%data
        deallocate(self%data)
        allocate(self%data(self%capacity * 2))
        self%data(1:self%capacity) = tmp
        self%capacity = self%capacity * 2
      end block
    end if
    self%top = self%top + 1
    self%data(self%top) = val
  end subroutine

  function stack_pop(self) result(val)
    class(stack_t), intent(inout) :: self
    integer :: val
    if (self%top == 0) stop "Stack underflow"
    val = self%data(self%top)
    self%top = self%top - 1
  end function

  function stack_peek(self) result(val)
    class(stack_t), intent(in) :: self
    integer :: val
    val = self%data(self%top)
  end function

  logical function stack_is_empty(self) result(r)
    class(stack_t), intent(in) :: self
    r = (self%top == 0)
  end function

  integer function stack_size(self) result(n)
    class(stack_t), intent(in) :: self
    n = self%top
  end function

  subroutine stack_clear(self)
    class(stack_t), intent(inout) :: self
    self%top = 0
  end subroutine

  subroutine stack_finalize(self)
    type(stack_t) :: self
    if (allocated(self%data)) deallocate(self%data)
  end subroutine
end module stack_mod

program use_stack
  use stack_mod
  implicit none
  type(stack_t) :: s
  s = stack_t(8)         ! initial capacity 8
  call s%push(10)
  call s%push(20)
  call s%push(30)
  print *, s%size()      ! 3
  print *, s%pop()       ! 30
  print *, s%pop()       ! 20
  print *, s%is_empty()  ! F
end program use_stack

module shapes
  implicit none
  private
  public :: shape, circle, rectangle, shape_ptr

  ! Base type (abstract)
  type, abstract :: shape
    character(20) :: name
  contains
    procedure(area_if), deferred :: area
    procedure(describe_if), deferred :: describe
    procedure :: get_name => shape_get_name
  end type shape

  ! Abstract interface (must be implemented by subclasses)
  abstract interface
    function area_if(self) result(a)
      import :: shape
      class(shape), intent(in) :: self
      real :: a
    end function area_if
    subroutine describe_if(self)
      import :: shape
      class(shape), intent(in) :: self
    end subroutine describe_if
  end interface

  ! Derived type: circle extends shape
  type, extends(shape) :: circle
    real :: radius
  contains
    procedure :: area => circle_area
    procedure :: describe => circle_describe
  end type circle

  ! Derived type: rectangle extends shape
  type, extends(shape) :: rectangle
    real :: width, height
  contains
    procedure :: area => rect_area
    procedure :: describe => rect_describe
  end type rectangle

  ! Class pointer for polymorphism
  type :: shape_ptr
    class(shape), allocatable :: ptr
  end type shape_ptr

contains
  function shape_get_name(self) result(n)
    class(shape), intent(in) :: self
    character(20) :: n
    n = self%name
  end function

  function circle_area(self) result(a)
    class(circle), intent(in) :: self
    real :: a
    a = 3.14159265 * self%radius**2
  end function

  subroutine circle_describe(self)
    class(circle), intent(in) :: self
    print *, "Circle '", trim(self%name), "' r=", self%radius
  end subroutine

  function rect_area(self) result(a)
    class(rectangle), intent(in) :: self
    real :: a
    a = self%width * self%height
  end function

  function rect_describe(self) result(s)
    class(rectangle), intent(in) :: self
    character(50) :: s
    write(s, '("Rectangle ", F6.2, "x", F6.2)') self%width, self%height
  end function
end module shapes

program use_polymorphism
  use shapes
  implicit none
  type(shape_ptr) :: shapes_arr(2)
  ! Allocate concrete types into polymorphic container
  allocate(circle::shapes_arr(1)%ptr)
  select type(s => shapes_arr(1)%ptr)
  type is (circle)
    s%radius = 5.0
    s%name = "C1"
  end select

  allocate(rectangle::shapes_arr(2)%ptr)
  select type(s => shapes_arr(2)%ptr)
  type is (rectangle)
    s%width = 3.0
    s%height = 4.0
    s%name = "R1"
  end select

  ! Polymorphic dispatch
  block
    integer :: i
    do i = 1, 2
      call shapes_arr(i)%ptr%describe
      print *, "Area: ", shapes_arr(i)%ptr%area()
    end do
  end block
end program use_polymorphism

program nested_types
  implicit none
  ! Nested derived types
  type :: Address
    character(50) :: street
    character(30) :: city
    character(10) :: zip
  end type Address

  type :: Employee
    integer :: id
    character(30) :: name
    type(Address) :: home      ! nested type
    real :: salary
  end type Employee

  type :: Department
    character(30) :: name
    type(Employee), allocatable :: employees(:)   ! array of types
    integer :: count = 0
  end type Department

  ! Construct with nested structure constructor
  type(Employee) :: emp
  emp = Employee(101, "Alice", &
        Address("123 Main St", "Springfield", "12345"), 75000.0)

  ! Access nested components
  print *, emp%name                      ! Alice
  print *, emp%home%city                 ! Springfield
  print *, emp%home%zip                  ! 12345

  ! Department with array of employees
  type(Department) :: dept
  dept%name = "Engineering"
  allocate(dept%employees(3))
  dept%employees(1) = emp
  dept%employees(2) = Employee(102, "Bob", &
                        Address("456 Oak Ave", "Shelbyville", "54321"), 68000.0)
  dept%count = 2

  ! Iterate over array of types
  block
    integer :: i
    do i = 1, dept%count
      print *, dept%employees(i)%id, trim(dept%employees(i)%name), &
               dept%employees(i)%home%city
    end do
  end block

  ! Array of derived type components
  print *, dept%employees(1:2)%id        ! array of id field
end program nested_types

program file_open
  implicit none
  integer :: u, ios
  character(100) :: msg

  ! newunit: compiler picks a unique unit number (Fortran 2008)
  ! status: 'old' (must exist), 'new' (must not exist),
  !         'replace' (overwrite), 'scratch' (temporary)
  ! action: 'read', 'write', 'readwrite'
  ! position: 'rewind', 'append', 'asis'
  open(newunit=u, file="data.txt", status="replace", &
       action="write", iostat=ios, iomsg=msg)
  if (ios /= 0) then
    print *, "Open failed: ", trim(msg)
    stop 1
  end if

  write(u, *) "First line"
  write(u, *) "Second line"
  close(u)

  ! Append to existing file
  open(newunit=u, file="data.txt", status="old", &
       action="write", position="append", iostat=ios)
  if (ios == 0) then
    write(u, *) "Appended line"
    close(u)
  end if

  ! Scratch file (auto-deleted on close)
  open(newunit=u, status="scratch", action="readwrite")
  write(u, *) "Temporary data"
  rewind(u)
  ! read back...
  close(u)   ! file disappears

  ! Reading with end-of-file detection
  open(newunit=u, file="data.txt", status="old", action="read")
  do
    read(u, '(A)', iostat=ios) msg
    if (ios /= 0) exit   ! EOF or error
    print *, trim(msg)
  end do
  close(u)
end program file_open

program formatted_io
  implicit none
  integer :: u, n = 42
  real :: pi = 3.14159265
  character(20) :: name = "Alice"

  ! Write formatted data to file
  open(newunit=u, file="output.txt", status="replace")
  write(u, '(A, I0)') "Count: ", n
  write(u, '(A, F8.4)') "Pi: ", pi
  write(u, '(A, A)') "Name: ", trim(name)
  write(u, '(3(I0, 1X))') 1, 2, 3   ! "1 2 3 "
  close(u)

  ! Read formatted data back
  open(newunit=u, file="output.txt", status="old", action="read")
  block
    character(100) :: line
    integer :: i
    do i = 1, 5
      read(u, '(A)') line
      print *, trim(line)
    end do
  end block
  close(u)

  ! Reading structured data
  open(newunit=u, file="data.csv", status="replace")
  write(u, '(I0, ",", I0, ",", I0)') 1, 10, 100
  write(u, '(I0, ",", I0, ",", I0)') 2, 20, 200
  close(u)

  ! Read back as numbers
  open(newunit=u, file="data.csv", status="old", action="read")
  block
    integer :: a, b, c
    ! List-directed read handles commas as separators
    read(u, *) a, b, c
    print *, a, b, c   ! 1 10 100
    read(u, *) a, b, c
    print *, a, b, c   ! 2 20 200
  end block
  close(u)
end program formatted_io

program binary_io
  implicit none
  integer :: u, i
  real :: arr(5) = [1.0, 2.0, 3.0, 4.0, 5.0]
  real :: arr_in(5)
  integer :: n

  ! Write binary (unformatted) - faster, no format conversion
  open(newunit=u, file="data.bin", status="replace", &
       form="unformatted", access="stream")
  write(u) size(arr)        ! write the count first
  write(u) arr              ! write the whole array
  close(u)

  ! Read binary back
  open(newunit=u, file="data.bin", status="old", &
       form="unformatted", access="stream", action="read")
  read(u) n                 ! read the count
  print *, "Count: ", n
  read(u) arr_in            ! read the array
  print *, arr_in           ! 1 2 3 4 5
  close(u)

  ! Sequential unformatted (default) - with record markers
  open(newunit=u, file="data_seq.bin", status="replace", &
       form="unformatted")   ! access defaults to "sequential"
  do i = 1, 5
    write(u) i, arr(i)      ! each write is a "record"
  end do
  close(u)

  ! Read sequentially
  open(newunit=u, file="data_seq.bin", status="old", &
       form="unformatted", action="read")
  do
    read(u, iostat=i) n, arr_in(1)
    if (i /= 0) exit
    print *, n, arr_in(1)
  end do
  close(u)
end program binary_io

program namelist_demo
  implicit none
  ! Namelist groups variables for easy I/O
  integer :: max_iter = 100
  real :: tolerance = 1.0e-6
  logical :: verbose = .true.
  character(20) :: method = "newton"
  real :: params(3) = [0.1, 0.2, 0.3]

  ! Define a namelist group
  namelist /config/ max_iter, tolerance, verbose, method, params

  ! Write namelist to file
  block
    integer :: u
    open(newunit=u, file="config.nml", status="replace")
    write(u, nml=config)
    close(u)
  end block

  ! The file looks like:
  ! &CONFIG
  !  MAX_ITER=100,
  !  TOLERANCE=1.0000000E-06,
  !  VERBOSE=T,
  !  METHOD="newton",
  !  PARAMS=0.100000, 0.200000, 0.300000,
  ! /

  ! Read namelist from file (overrides defaults)
  block
    integer :: u, ios
    open(newunit=u, file="config.nml", status="old", action="read")
    read(u, nml=config, iostat=ios)
    close(u)
    if (ios /= 0) then
      print *, "Error reading config"
    else
      print *, "max_iter=", max_iter
      print *, "tolerance=", tolerance
      print *, "method=", trim(method)
      print *, "params=", params
    end if
  end block

  ! User can edit config.nml in a text editor
  ! then re-run to pick up new values
end program namelist_demo

program internal_io
  implicit none
  character(100) :: buffer
  integer :: n = 42
  real :: x = 3.14159
  integer :: ios

  ! Internal write: format to a string (like sprintf)
  write(buffer, '(A, I0, A, F8.4)') "n=", n, " x=", x
  print *, trim(buffer)   ! n=42 x=  3.1416

  ! Internal read: parse from a string (like sscanf)
  character(50) :: input = "100 3.14 hello"
  integer :: a
  real :: b
  character(20) :: c
  read(input, *, iostat=ios) a, b, c
  if (ios == 0) print *, a, b, trim(c)   ! 100 3.14 hello

  ! Robust number parsing with error handling
  character(20) :: num_str = "3.14abc"
  real :: val
  read(num_str, *, iostat=ios) val
  if (ios /= 0) then
    print *, "Parse error: '", trim(num_str), "' is not a number"
  else
    print *, "Value: ", val
  end if

  ! End-of-file and error handling on file reads
  block
    integer :: u
    character(100) :: line
    open(newunit=u, file="data.txt", status="old", action="read")
    do
      read(u, '(A)', iostat=ios) line
      if (ios < 0) then
        print *, "End of file"
        exit
      else if (ios > 0) then
        print *, "Read error at line"
        exit
      end if
      print *, trim(line)
    end do
    close(u)
  end block

  ! err= and end= labels (legacy style)
  block
    integer :: u, val
    open(newunit=u, file="nums.txt", status="old", action="read")
    do
      read(u, *, end=100, err=200) val
      print *, val
    end do
100 print *, "Reached EOF"
    close(u)
    goto 300
200 print *, "Read error!"
    close(u)
300 continue
  end block
end program internal_io

program precision_demo
  use iso_fortran_env, only: int32, int64, real32, real64, real128
  implicit none
  ! Portable kind selection via iso_fortran_env
  integer(kind=int32) :: i32 = 100
  integer(kind=int64) :: i64 = 9223372036854775807_int64
  real(kind=real32) :: r32 = 3.14159265_real32   ! single (~7 digits)
  real(kind=real64) :: r64 = 3.14159265358979_real64  ! double (~15 digits)
  real(kind=real128) :: r128 = 3.14159265358979_real128 ! quad (~33 digits)

  ! Legacy: selected_real_kind (still works)
  integer, parameter :: dp = selected_real_kind(15, 307)  ! double
  integer, parameter :: sp = selected_real_kind(6, 37)    ! single
  real(kind=dp) :: pi = 3.14159265358979_dp

  print *, "real32 precision: ", precision(r32), " range: ", range(r32)
  print *, "real64 precision: ", precision(r64), " range: ", range(r64)
  print *, "real128 precision: ", precision(r128)

  ! Epsilon and tiny
  print *, "epsilon(r32): ", epsilon(r32)   ! ~1.19e-7
  print *, "epsilon(r64): ", epsilon(r64)   ! ~2.22e-16
  print *, "tiny(r64): ", tiny(r64)         ! smallest positive
  print *, "huge(i32): ", huge(i32)         ! 2147483647
  print *, "huge(i64): ", huge(i64)

  ! Always use _kind suffix on literals
  ! WRONG: real(kind=dp) :: x = 3.14  (loses precision!)
  ! RIGHT: real(kind=dp) :: x = 3.14_dp
end program precision_demo

program linalg
  implicit none
  real(8) :: A(3,3), B(3,3), C(3,3)
  real(8) :: v(3), w(3), x(3)
  real(8) :: det
  integer :: i

  ! Initialize matrices
  A = reshape([1,2,3, 4,5,6, 7,8,10], [3,3])  ! column-major fill
  B = reshape([1,0,0, 0,1,0, 0,0,1], [3,3])   ! identity

  ! Matrix-matrix multiplication
  C = matmul(A, B)        ! A * I = A
  print *, "A * I = A? ", all(abs(C - A) < 1e-10)

  ! Matrix-vector multiplication
  v = [1.0, 2.0, 3.0]
  w = matmul(A, v)
  print *, "A * v = ", w    ! 14 32 53

  ! Dot product
  print *, "v . v = ", dot_product(v, v)   ! 14

  ! Outer product
  C = 0.0
  do i = 1, 3
    C(:,i) = v * v(i)
  end do

  ! Transpose
  C = transpose(A)

  ! Solve linear system Ax = b (need LAPACK or custom)
  ! Using LAPACK dgesv:
  !   call dgesv(n, nrhs, A, lda, ipiv, b, ldb, info)
  ! For demo, manual Gaussian elimination:
  x = solve_3x3(A, v)
  print *, "Solution: ", x

contains
  function solve_3x3(A, b) result(x)
    real(8), intent(inout) :: A(3,3)
    real(8), intent(in) :: b(3)
    real(8) :: x(3), M(3,4), factor
    integer :: k, j
    M(:,1:3) = A
    M(:,4) = b
    ! Forward elimination
    do k = 1, 2
      do i = k+1, 3
        factor = M(i,k) / M(k,k)
        M(i,:) = M(i,:) - factor * M(k,:)
      end do
    end do
    ! Back substitution
    x(3) = M(3,4) / M(3,3)
    x(2) = (M(2,4) - M(2,3)*x(3)) / M(2,2)
    x(1) = (M(1,4) - M(1,2)*x(2) - M(1,3)*x(3)) / M(1,1)
  end function
end program linalg

program random_demo
  implicit none
  real :: r
  integer :: i, n = 10
  real :: arr(10)
  integer :: seed_size
  integer, allocatable :: seed(:)

  ! Initialize random seed (Fortran 2003: random_seed with no args
  ! uses a processor-dependent seed; for reproducibility, set it)
  call random_seed(size=seed_size)
  allocate(seed(seed_size))
  ! Set a fixed seed for reproducibility
  seed = [(12345 + i*6789, i=1, seed_size)]
  call random_seed(put=seed)

  ! Generate uniform [0,1) random reals
  call random_number(r)
  print *, "Single: ", r

  call random_number(arr)   ! fills entire array
  print *, "Array: ", arr

  ! Generate integers in [a, b]
  block
    integer :: a = 1, b = 6
    integer :: dice
    do i = 1, 5
      call random_number(r)
      dice = a + int(r * (b - a + 1))   ! [1, 6]
      print *, "Dice roll: ", dice
    end do
  end block

  ! Normal (Gaussian) via Box-Muller transform
  block
    real :: u1, u2, z1, z2
    integer :: j
    do j = 1, 5
      call random_number(u1)
      call random_number(u2)
      u1 = max(u1, 1e-10)   ! avoid log(0)
      z1 = sqrt(-2.0 * log(u1)) * cos(2.0 * 3.14159265 * u2)
      z2 = sqrt(-2.0 * log(u1)) * sin(2.0 * 3.14159265 * u2)
      print *, "Gaussian: ", z1, z2
    end do
  end block

  ! Monte Carlo: estimate pi
  block
    integer :: inside = 0, total = 100000
    real :: x, y
    do i = 1, total
      call random_number(x)
      call random_number(y)
      if (x*x + y*y <= 1.0) inside = inside + 1
    end do
    print *, "Pi estimate: ", 4.0 * real(inside) / total
  end block
end program random_demo

program numerical
  implicit none
  real(8) :: a, b, result
  integer :: n

  ! Numerical integration: Simpson's rule
  ! Integrate f(x) = x^2 from 0 to 2 (exact: 8/3 ≈ 2.6667)
  a = 0.0; b = 2.0; n = 1000
  result = simpson(f_sq, a, b, n)
  print *, "Integral of x^2 from 0 to 2: ", result   ! ~2.6667

  ! Integrate sin(x) from 0 to pi (exact: 2.0)
  result = simpson(f_sin, 0.0_8, 3.14159265358979_8, 1000)
  print *, "Integral of sin(x) from 0 to pi: ", result  ! ~2.0

  ! Root finding: bisection method
  ! Find root of f(x) = x^2 - 2 (i.e., sqrt(2) ≈ 1.4142)
  result = bisection(f_x2_minus_2, 0.0_8, 2.0_8, 1e-12_8)
  print *, "sqrt(2) = ", result   ! ~1.41421356

  ! Newton's method (needs derivative)
  result = newton(f_x2_minus_2, fp_x2_minus_2, 1.0_8, 1e-12_8)
  print *, "sqrt(2) via Newton: ", result

contains
  ! Function to integrate: x^2
  function f_sq(x) result(y)
    real(8), intent(in) :: x
    real(8) :: y
    y = x * x
  end function

  function f_sin(x) result(y)
    real(8), intent(in) :: x
    real(8) :: y
    y = sin(x)
  end function

  ! Simpson's rule: integral of f from a to b with n intervals
  function simpson(f, a, b, n) result(integral)
    interface
      function f(x) result(y)
        import
        real(8), intent(in) :: x
        real(8) :: y
      end function
    end interface
    real(8), intent(in) :: a, b
    integer, intent(in) :: n
    real(8) :: integral, h, x
    integer :: i
    h = (b - a) / n
    integral = f(a) + f(b)
    do i = 1, n-1
      x = a + i * h
      if (mod(i, 2) == 0) then
        integral = integral + 2.0 * f(x)
      else
        integral = integral + 4.0 * f(x)
      end if
    end do
    integral = integral * h / 3.0
  end function

  function f_x2_minus_2(x) result(y)
    real(8), intent(in) :: x
    real(8) :: y
    y = x*x - 2.0
  end function

  function fp_x2_minus_2(x) result(y)
    real(8), intent(in) :: x
    real(8) :: y
    y = 2.0 * x
  end function

  ! Bisection: find root in [a, b] (f(a) and f(b) must have opposite signs)
  function bisection(f, a, b, tol) result(root)
    interface
      function f(x) result(y)
        import
        real(8), intent(in) :: x
        real(8) :: y
      end function
    end interface
    real(8), intent(in) :: a, b, tol
    real(8) :: root, fa, fb, mid, fmid
    integer :: iter, maxiter = 100
    fa = f(a); fb = f(b)
    if (fa * fb > 0) stop "No sign change in interval"
    do iter = 1, maxiter
      mid = (a + b) / 2.0
      fmid = f(mid)
      if (abs(fmid) < tol .or. (b - a)/2.0 < tol) then
        root = mid
        return
      end if
      if (fa * fmid < 0) then
        b = mid; fb = fmid
      else
        a = mid; fa = fmid
      end if
    end do
    root = (a + b) / 2.0
  end function

  ! Newton's method: x_{n+1} = x_n - f(x)/f'(x)
  function newton(f, fp, x0, tol) result(root)
    interface
      function f(x) result(y)
        import
        real(8), intent(in) :: x
        real(8) :: y
      end function
      function fp(x) result(y)
        import
        real(8), intent(in) :: x
        real(8) :: y
      end function
    end interface
    real(8), intent(in) :: x0, tol
    real(8) :: root, x, fx, fpx
    integer :: iter, maxiter = 100
    x = x0
    do iter = 1, maxiter
      fx = f(x); fpx = fp(x)
      if (abs(fx) < tol) then
        root = x
        return
      end if
      x = x - fx / fpx
    end do
    root = x
  end function
end program numerical

program ieee_demo
  use ieee_arithmetic
  implicit none
  real :: a, b, c
  logical :: flag

  ! IEEE special values
  a = ieee_value(a, ieee_positive_inf)   ! +Infinity
  b = ieee_value(b, ieee_negative_inf)   ! -Infinity
  c = ieee_value(c, ieee_quiet_nan)      ! NaN

  print *, "Infinity: ", a                ! Infinity
  print *, "NaN: ", c                     ! NaN
  print *, "Inf > 1e30: ", a > 1e30       ! T
  print *, "NaN == NaN: ", c == c         ! F (NaN is never equal!)

  ! Check for special values
  print *, "is_nan(c): ", ieee_is_nan(c)         ! T
  print *, "is_finite(1.0): ", ieee_is_finite(1.0)  ! T
  print *, "is_inf(a): ", ieee_is_finite(a)      ! F

  ! Operations producing special values
  print *, "1.0/0.0: ", 1.0/0.0           ! Infinity (if -ffast-math off)
  print *, "0.0/0.0: ", 0.0/0.0           ! NaN
  print *, "sqrt(-1.0): ", sqrt(-1.0)     ! NaN

  ! IEEE exception flags
  call ieee_set_halting_mode(ieee_divide_by_zero, .false.)  ! don't halt
  b = 1.0 / 0.0   ! sets divide_by_zero flag, returns Inf
  call ieee_get_flag(ieee_divide_by_zero, flag)
  print *, "Divide by zero occurred: ", flag   ! T

  ! Check and clear flags
  call ieee_set_flag(ieee_all, .false.)   ! clear all flags
  b = 1.0 / 0.0
  call ieee_get_flag(ieee_divide_by_zero, flag)
  print *, "Flag after division: ", flag   ! T
  call ieee_set_flag(ieee_all, .false.)   ! clear

  ! Rounding modes
  call ieee_set_rounding_mode(ieee_nearest)   ! default
  call ieee_set_rounding_mode(ieee_down)      ! round toward -inf
  call ieee_set_rounding_mode(ieee_up)        ! round toward +inf
  call ieee_set_rounding_mode(ieee_to_zero)   ! truncate

  ! Comparing NaN-safe
  if (ieee_unordered(c, 1.0)) print *, "c is unordered (NaN)"
end program ieee_demo

program coarray_hello
  use iso_fortran_env, only: real64
  implicit none
  real(real64) :: x[*]  ! coarray — one copy per image
  integer :: me

  me = this_image()
  x = real(me, real64)  ! local assignment

  call co_sum(x, result_image=1)  ! reduce sum to image 1
  if (me == 1) print *, 'Sum =', x

  sync all  ! barrier
end program

program remote_access
  implicit none
  integer :: val[*], neighbor
  integer :: me, n

  me = this_image()
  n = num_images()
  val = me * 10

  sync all  ! ensure all writes complete

  ! read from neighbor (circular)
  neighbor = merge(1, me + 1, me == n)
  print *, 'Image', me, 'sees neighbor', neighbor, 'value', val[neighbor]

  ! write to image 1 from all
  if (me /= 1) val[1] = val[me]
  sync all
  if (me == 1) print *, 'Image 1 received:', val
end program

program collectives
  use iso_fortran_env, only: real64
  implicit none
  real(real64) :: local_sum, global_sum[*]
  integer :: i

  local_sum = 0.0_real64
  do i = 1, 100
    local_sum = local_sum + real(i * this_image(), real64)
  end do

  global_sum = local_sum
  call co_sum(global_sum, result_image=1)
  if (this_image() == 1) print *, 'Total:', global_sum

  ! other collectives: co_min, co_max, co_broadcast
  call co_broadcast(global_sum, source_image=1)
end program

program coarray_types
  implicit none
  type :: particle
    real :: x, y, z
    real :: mass
  end type
  type(particle) :: p[*]
  integer :: me

  me = this_image()
  p%mass = real(me)
  p%x = real(me) * 0.5

  sync all
  ! access component of remote coarray
  if (me == 1) print *, 'Image 2 mass:', p[2]%mass
end program

program alloc_coarrays
  use iso_fortran_env, only: team_type
  implicit none
  integer, allocatable :: data(:)[:]
  type(team_type) :: odd_team, even_team
  integer :: me

  me = this_image()
  allocate(data(100)[*])  ! allocate on all images

  ! F2018 teams — split images into groups
  form team(merge(1, 2, mod(me, 2) == 1), odd_team, even_team)
  if (mod(me, 2) == 1) then
    change team(odd_team)
      ! this_image() and num_images() now refer to team
      data(this_image()) = me
    end team
  end if

  deallocate(data)
end program

program c_interop
  use iso_c_binding, only: c_int, c_double, c_char, c_null_char
  implicit none
  interface
    subroutine c_print(msg) bind(C, name='c_print')
      import :: c_char
      character(kind=c_char), dimension(*) :: msg
    end subroutine
  end interface

  call c_print('Hello from Fortran' // c_null_char)
end program

! Corresponding C:
!   void c_print(const char* msg) { printf("%s\n", msg); }

program array_interop
  use iso_c_binding, only: c_double, c_int, c_loc, c_f_pointer
  implicit none
  real(c_double), target :: arr(10)
  integer(c_int) :: n
  type(c_ptr) :: ptr

  arr = [(real(i), i=1,10)]
  ptr = c_loc(arr(1))  ! get C pointer

  ! call C function: void process(double* arr, int n);
  call process_c(ptr, size(arr, kind=c_int))

  interface
    subroutine process_c(arr, n) bind(C, name='process')
      import :: c_double, c_int, c_ptr
      type(c_ptr), value :: arr
      integer(c_int), value :: n
    end subroutine
  end interface
end program

module data_types
  use iso_c_binding, only: c_double, c_int, c_char
  implicit none

  type, bind(C) :: point
    real(c_double) :: x, y
    integer(c_int) :: id
  end type

  type, bind(C) :: string_holder
    character(kind=c_char, len=1) :: name(64)
  end type
end module

! C equivalent:
!   struct point { double x, y; int id; };
!   struct string_holder { char name[64]; };

! Fortran:
module fmod
  use iso_c_binding, only: c_double
  implicit none
contains
  subroutine square_array(arr, n) bind(C, name='square_array')
    integer, value :: n
    real(c_double), intent(inout) :: arr(n)
    arr = arr ** 2
  end subroutine
end module

! C caller:
!   extern void square_array(double* arr, int n);
!   double data[5] = {1, 2, 3, 4, 5};
!   square_array(data, 5);

module callbacks
  use iso_c_binding, only: c_funloc, c_funptr, c_int
  implicit none

  abstract interface
    function comparator(a, b) bind(C)
      import :: c_int
      integer(c_int), value :: a, b
      integer(c_int) :: comparator
    end function
  end interface

contains
  subroutine sort_with_c(arr, n, cmp) bind(C)
    integer(c_int), value :: n
    integer(c_int), intent(inout) :: arr(n)
    procedure(comparator) :: cmp
    ! ... use cmp(a, b) to compare ...
  end subroutine
end module

! C side:
!   int descending(int a, int b) { return b - a; }
!   sort_with_c(arr, n, descending);

module shapes
  implicit none
  type, abstract :: shape
    real :: x = 0, y = 0
  contains
    procedure(area_iface), deferred :: area
    procedure :: move => shape_move
  end type

  abstract interface
    function area_iface(this) result(a)
      import :: shape
      class(shape), intent(in) :: this
      real :: a
    end function
  end interface

  type, extends(shape) :: circle
    real :: radius
  contains
    procedure :: area => circle_area
  end type

contains
  subroutine shape_move(this, dx, dy)
    class(shape), intent(inout) :: this
    real, intent(in) :: dx, dy
    this%x = this%x + dx
    this%y = this%y + dy
  end subroutine

  function circle_area(this) result(a)
    class(circle), intent(in) :: this
    real :: a
    a = 3.14159 * this%radius ** 2
  end function
end module

module polymorph
  use shapes, only: shape, circle, rectangle
  implicit none
contains
  subroutine describe(s)
    class(shape), intent(in) :: s
    select type(s)
    type is (circle)
      print *, 'Circle radius:', s%radius
    type is (rectangle)
      print *, 'Rectangle:', s%width, 'x', s%height
    class is (shape)
      print *, 'Some shape at:', s%x, s%y
    class default
      print *, 'Unknown type'
    end select
  end subroutine

  subroutine process_all(shapes)
    class(shape), intent(in) :: shapes(:)
    integer :: i
    do i = 1, size(shapes)
      call describe(shapes(i))
      print *, 'Area:', shapes(i)%area()
    end do
  end subroutine
end module

module resources
  implicit none
  type :: file_handle
    integer :: unit = -1
  contains
    final :: close_file
    procedure :: open
  end type
contains
  subroutine open(this, filename)
    class(file_handle), intent(inout) :: this
    character(*), intent(in) :: filename
    open(newunit=this%unit, file=filename, status='old')
  end subroutine

  subroutine close_file(this)
    type(file_handle), intent(inout) :: this
    if (this%unit /= -1) then
      close(this%unit)
      this%unit = -1
    end if
  end subroutine
end module

program demo
  use resources
  type(file_handle) :: f
  call f%open('data.txt')
  ! ... use file ...
end program  ! f goes out of scope -> close_file called automatically

module points
  implicit none
  type :: point
    real :: x, y
  contains
    procedure :: norm
  end type
  interface point
    procedure new_point
    procedure new_point_polar
  end interface
contains
  function new_point(x, y) result(p)
    real, intent(in) :: x, y
    type(point) :: p
    p%x = x
    p%y = y
  end function

  function new_point_polar(r, theta) result(p)
    real, intent(in) :: r, theta
    type(point) :: p
    p%x = r * cos(theta)
    p%y = r * sin(theta)
  end function

  function norm(this) result(n)
    class(point), intent(in) :: this
    real :: n
    n = sqrt(this%x**2 + this%y**2)
  end function
end module

program use_points
  use points
  type(point) :: a, b
  a = point(1.0, 2.0)         ! cartesian
  b = point(3.0, 0.5)         ! polar (same name, different args)
  print *, a%norm(), b%norm()
end program

module container
  implicit none
  type, abstract :: container
  contains
    procedure(add_iface), deferred :: add
    procedure(size_iface), deferred :: size
    procedure :: is_empty => container_is_empty
  end type

  abstract interface
    subroutine add_iface(this, item)
      import :: container
      class(container), intent(inout) :: this
      class(*), intent(in) :: item
    end subroutine
    function size_iface(this) result(n)
      import :: container
      class(container), intent(in) :: this
      integer :: n
    end function
  end interface
contains
  function container_is_empty(this) result(e)
    class(container), intent(in) :: this
    logical :: e
    e = (this%size() == 0)
  end function
end module

! Concrete subclass implements add/size
! class(*) allows storing any type (unlimited polymorphic)

module pdt_types
  implicit none
  type :: matrix(k, n, m)
    integer, kind :: k = kind(1.0)
    integer, len :: n, m
    real(k) :: data(n, m)
  end type

  ! kind parameters: compile-time (like template params)
  ! len parameters: runtime (like allocatable dimensions)
end module

program use_pdt
  use pdt_types
  implicit none
  type(matrix(kind(1.0), 3, 3)) :: small  ! fixed 3x3
  type(matrix(kind(1.d0), :, :)), allocatable :: big  ! deferred

  small%data = 0.0
  allocate(matrix(kind(1.d0), 100, 100) :: big)
  big%data = 0.d0
  deallocate(big)
end program

module typed_array
  implicit none
  type :: array_t(k)
    integer, kind :: k = kind(1.0)
    real(k), allocatable :: data(:)
  contains
    procedure :: sum_values
  end type
contains
  function sum_values(this) result(s)
    class(array_t(k=*)), intent(in) :: this
    real(this%k) :: s
    s = sum(this%data)
  end function
end module

program test
  use typed_array
  type(array_t(kind(1.0)))  :: float_arr
  type(array_t(kind(1.d0))) :: double_arr

  allocate(float_arr%data(10))
  allocate(double_arr%data(10))
  float_arr%data = 1.0
  double_arr%data = 1.d0

  print *, float_arr%sum_values()
  print *, double_arr%sum_values()
end program

module fixed_string
  implicit none
  type :: fstr(n)
    integer, len :: n
    character(len=n) :: str
  end type
contains
  function concat(a, b) result(c)
    type(fstr(*)), intent(in) :: a, b
    type(fstr(a%n + b%n)) :: c
    c%str = a%str // b%str
  end function
end module

program test
  use fixed_string
  type(fstr(5))  :: a = fstr(5)('Hello')
  type(fstr(6))  :: b = fstr(6)(' World')
  type(fstr(11)) :: c

  c = concat(a, b)
  print *, c%str
end program

module stack_t
  implicit none
  type :: stack(k)
    integer, kind :: k = kind(1.0)
    real(k), allocatable :: data(:)
    integer :: top = 0
  contains
    procedure :: push
    procedure :: pop
    procedure :: is_empty
  end type
contains
  subroutine push(this, val)
    class(stack(k=*)), intent(inout) :: this
    real(this%k), intent(in) :: val
    if (this%top == size(this%data)) then
      block
        real(this%k), allocatable :: tmp(:)
        tmp = [this%data, val]
        call move_alloc(tmp, this%data)
      end block
    else
      this%top = this%top + 1
      this%data(this%top) = val
    end if
  end subroutine
end module

! Limitations of PDTs:
! - Not all compilers fully support len parameters
! - I/O of PDTs may not work as expected
! - Some F2003 features (e.g., PDT arrays) have spotty support

! Workaround: use allocatable components instead of len params
module alt_matrix
  implicit none
  type :: matrix_alt(k)
    integer, kind :: k = kind(1.0)
    real(k), allocatable :: data(:,:)
  end type
contains
  function make_matrix(k_, n, m) result(mat)
    integer, intent(in) :: n, m
    ! kind must be compile-time constant — can't pass dynamically
    ! workaround: separate constructors per kind
  end function
end module

! Best practice: use kind params (well-supported), avoid len params
! Use allocatable components for runtime sizing instead

! parent_module.f90
module math_lib
  implicit none
  interface
    module function integrate(f, a, b, n) result(s)
      abstract interface
        real function func_t(x)
          real, intent(in) :: x
        end function
      end interface
      procedure(func_t) :: f
      real, intent(in) :: a, b
      integer, intent(in) :: n
      real :: s
    end function
  end interface
end module

! sub_math.f90
submodule(math_lib) math_impl
  implicit none
contains
  module procedure integrate
    integer :: i
    real :: dx, x
    dx = (b - a) / n
    s = 0.5 * (f(a) + f(b))
    do i = 1, n-1
      x = a + i * dx
      s = s + f(x)
    end do
    s = s * dx
  end procedure
end submodule

! big_lib.f90
module big_lib
  implicit none
  interface
    module subroutine sort(arr)
      real, intent(inout) :: arr(:)
    end subroutine
    module function mean(arr) result(m)
      real, intent(in) :: arr(:)
      real :: m
    end function
    module function stddev(arr) result(s)
      real, intent(in) :: arr(:)
      real :: s
    end function
  end interface
end module

! sub_sort.f90 — one submodule
submodule(big_lib) sort_impl
contains
  module procedure sort
    ! quicksort implementation
  end procedure
end submodule

! sub_stats.f90 — another submodule
submodule(big_lib) stats_impl
contains
  module procedure mean
    m = sum(arr) / size(arr)
  end procedure
  module procedure stddev
    block
      real :: mu
      mu = mean(arr)  ! can call sibling module procedures
      s = sqrt(sum((arr - mu)**2) / (size(arr) - 1))
    end block
  end procedure
end submodule

module sparse_ops
  implicit none
  interface
    module function sparse_multiply(a, b) result(c)
      ! ... sparse matrix types ...
      real, allocatable :: c(:,:)
    end function
  end interface
end module

submodule(sparse_ops) sparse_impl
  implicit none
contains
  module procedure sparse_multiply
    ! call helper — invisible to module users
    call validate_dimensions(a, b)
    c = multiply_kernel(a, b)
  end procedure

  ! private helper — not exposed via parent module
  subroutine validate_dimensions(a, b)
    ! ...
  end subroutine

  function multiply_kernel(a, b) result(c)
    ! ...
  end function
end submodule

module config
  implicit none
  integer :: verbosity = 0
  interface
    module subroutine log_message(msg)
      character(*), intent(in) :: msg
    end subroutine
  end interface
end module

submodule(config) config_impl
  implicit none
  ! submodule can access parent module's variables
contains
  module procedure log_message
    if (verbosity > 0) then
      print *, '[LOG]', msg
    end if
  end procedure
end submodule

! Usage:
!   use config
!   verbosity = 1
!   call log_message('Starting up')

! Without submodules:
!   module big_mod
!     ... 5000 lines of implementation ...
!   end module
!   ! Any change -> recompile big_mod + all users

! With submodules:
!   module big_mod
!     ! just interfaces (~200 lines)
!   end module
!   submodule(big_mod) impl
!     ! 5000 lines of implementation
!   end submodule
!   ! Change impl body -> recompile only submodule
!   ! Change interface -> recompile everything (unavoidable)

! Build system example (Makefile):
!   big_mod.o: big_mod.f90
!   impl.o: impl.f90 big_mod.o
!   user.o: user.f90 big_mod.o  ! NOT impl.o
!   app: user.o impl.o big_mod.o
!   	$(FC) -o app user.o impl.o big_mod.o

program ieee_exceptions
  use ieee_exceptions
  use ieee_arithmetic, only: ieee_is_nan, ieee_is_finite
  implicit none
  real :: x, y
  logical :: overflow_flag

  call ieee_set_halting_mode(ieee_overflow, .false.)
  x = huge(x) * 10  ! overflow -> Inf, no halt

  call ieee_get_flag(ieee_overflow, overflow_flag)
  print *, 'Overflow occurred:', overflow_flag

  y = 0.0 / 0.0  ! NaN
  print *, 'x is finite:', ieee_is_finite(x)
  print *, 'y is NaN:', ieee_is_nan(y)

  call ieee_set_flag(ieee_all, .false.)  ! clear all flags
end program

program nan_inf
  use ieee_arithmetic, only: ieee_value, ieee_nan, &
                             ieee_positive_inf, ieee_negative_inf
  implicit none
  real :: pos_inf, neg_inf, nan_val

  pos_inf = ieee_value(0.0, ieee_positive_inf)
  neg_inf = ieee_value(0.0, ieee_negative_inf)
  nan_val = ieee_value(0.0, ieee_nan)

  print *, 'Inf - Inf =', pos_inf - pos_inf  ! NaN
  print *, '0 * Inf =', 0.0 * pos_inf         ! NaN
  print *, 'Inf > 1e30:', pos_inf > 1e30       ! T
  print *, 'NaN == NaN:', nan_val == nan_val   ! F (always!)

  ! generate via arithmetic
  pos_inf = 1.0 / 0.0  ! may halt without -fno-trapping-math
end program

program rounding_modes
  use ieee_arithmetic, only: ieee_set_rounding_mode, &
                             ieee_nearest, ieee_down, ieee_up, ieee_to_zero
  implicit none
  real :: x, y, result

  x = 1.0 / 3.0
  print *, 'Nearest:', x

  call ieee_set_rounding_mode(ieee_down)
  result = 1.0 / 3.0
  print *, 'Down:', result

  call ieee_set_rounding_mode(ieee_up)
  result = 1.0 / 3.0
  print *, 'Up:', result

  call ieee_set_rounding_mode(ieee_to_zero)
  result = 1.0 / 3.0
  print *, 'To zero:', result

  call ieee_set_rounding_mode(ieee_nearest)  ! restore default
end program

program ieee_inquiry
  use ieee_features
  use ieee_arithmetic, only: ieee_support_nan, ieee_support_inf, &
                             ieee_support_rounding, ieee_nearest
  implicit none

  if (ieee_support_nan(1.0)) print *, 'NaN supported'
  if (ieee_support_inf(1.0)) print *, 'Inf supported'
  if (ieee_support_rounding(ieee_nearest, 1.0)) &
    print *, 'Nearest rounding supported'

  ! Check halting mode
  block
    use ieee_exceptions
    logical :: halting
    call ieee_get_halting_mode(ieee_overflow, halting)
    print *, 'Overflow halts:', halting
  end block
end program

program fp_env
  use ieee_arithmetic
  use ieee_exceptions
  implicit none
  real :: a, b, c

  ! Save/restore FP state
  block
    use iso_fortran_env, only: real32
    real(real32) :: saved_state
    ! (no portable save/restore in standard — use compiler intrinsics)
  end block

  ! Disable halting for a risky section
  call ieee_set_halting_mode(ieee_all, .false.)
  call ieee_set_flag(ieee_all, .false.)

  a = sqrt(-1.0)  ! NaN, no halt
  b = 1.0 / 0.0   ! Inf, no halt

  if (ieee_is_nan(a)) then
    print *, 'Got NaN, using fallback'
    a = 0.0
  end if

  ! Re-enable halting
  call ieee_set_halting_mode(ieee_divide_by_zero, .true.)
  call ieee_set_halting_mode(ieee_invalid, .true.)
end program

program array_perf
  implicit none
  real, allocatable :: a(:,:), b(:,:)
  integer :: i, j, n
  real :: s

  n = 1000
  allocate(a(n,n), b(n,n))

  ! GOOD: column-major access (Fortran is column-major)
  do j = 1, n
    do i = 1, n
      s = s + a(i,j)
    end do
  end do

  ! BAD: row-major access (cache misses)
  do i = 1, n
    do j = 1, n
      s = s + a(i,j)
    end do
  end do

  ! BEST: array syntax (compiler optimizes)
  s = sum(a)
end program

module math_ops
  implicit none
contains
  ! pure: no side effects, enables optimization
  pure function square(x) result(y)
    real, intent(in) :: x
    real :: y
    y = x * x
  end function

  ! elemental: works on scalars and arrays element-wise
  pure elemental function clamp(x, lo, hi) result(y)
    real, intent(in) :: x, lo, hi
    real :: y
    y = max(lo, min(x, hi))
  end function
end module

program use_math
  use math_ops
  implicit none
  real :: arr(100), scalar

  arr = clamp(arr, 0.0, 1.0)  ! applies to each element
  scalar = clamp(1.5, 0.0, 1.0)  ! also works on scalar
  arr = square(arr)  ! elemental too (if declared)
end program

! gfortran optimization flags:
!   -O0   no optimization (debug)
!   -O1   basic optimization
!   -O2   standard optimization (recommended)
!   -O3   aggressive (may increase code size)
!   -Ofast  -O3 + -ffast-math (breaks IEEE compliance)
!   -march=native  use CPU's full instruction set
!   -flto  link-time optimization (cross-file inlining)
!   -fopenmp  enable OpenMP pragmas
!   -fprofile-generate/use  profile-guided optimization

! Example Makefile:
!   FC = gfortran
!   FFLAGS = -O2 -march=native -Wall -fcheck=all
!   FFLAGS_RELEASE = -O3 -march=native -flto -fno-trapping-math
!   FFLAGS_DEBUG = -O0 -g -fcheck=all -fbacktrace -ffpe-trap=invalid,zero

! ifort equivalents:
!   -O2, -O3, -xHost (= -march=native), -ipo (= -flto)

program openmp_demo
  use omp_lib
  implicit none
  real, allocatable :: a(:), b(:), c(:)
  integer :: i, n, nthreads

  n = 10000000
  allocate(a(n), b(n), c(n))
  a = 1.0; b = 2.0

  !$omp parallel
  if (omp_get_thread_num() == 0) then
    nthreads = omp_get_num_threads()
    print *, 'Threads:', nthreads
  end if
  !$omp end parallel

  !$omp parallel do
  do i = 1, n
    c(i) = a(i) + b(i)
  end do
  !$omp end parallel do

  ! reduction
  block
    real :: s
    s = 0.0
    !$omp parallel do reduction(+:s)
    do i = 1, n
      s = s + c(i)
    end do
    !$omp end parallel do
    print *, 'Sum:', s
  end block
end program

! Compile with profiling:
!   gfortran -pg -O2 prog.f90 -o prog
!   ./prog  (generates gmon.out)
!   gprof prog gmon.out > profile.txt

! Or use perf (Linux):
!   perf record ./prog
!   perf report

! Or gprofng (modern):
!   gprofng collect app ./prog
!   gprofng analyze

! Manual timing:
program timing
  use iso_fortran_env, only: real64, int64
  implicit none
  integer(int64) :: start, finish, rate
  real(real64) :: elapsed
  real, allocatable :: a(:)
  integer :: i

  allocate(a(10000000))
  call system_clock(start, rate)
  do i = 1, 100
    a = a + 1.0
  end do
  call system_clock(finish)
  elapsed = real(finish - start, real64) / real(rate, real64)
  print *, 'Elapsed:', elapsed, 'sec'
end program

module c_glue
  use iso_c_binding
  implicit none
  interface
    function c_malloc(size) bind(C, name='malloc')
      import :: c_ptr, c_size_t
      integer(c_size_t), value :: size
      type(c_ptr) :: c_malloc
    end function

    subroutine c_free(ptr) bind(C, name='free')
      import :: c_ptr
      type(c_ptr), value :: ptr
    end subroutine

    subroutine c_qsort(base, nmemb, size, compar) bind(C, name='qsort')
      import :: c_ptr, c_size_t, c_int
      type(c_ptr), value :: base
      integer(c_size_t), value :: nmemb, size
      abstract interface
        function compar_fn(a, b) bind(C)
          import :: c_ptr
          type(c_ptr), value :: a, b
          integer(c_int) :: compar_fn
        end function
      end interface
      procedure(compar_fn) :: compar
    end subroutine
  end interface
end module

! C++ side (wrapper.cpp):
!   extern "C" {
!     void cpp_process(double* data, int n) {
!       std::vector<double> v(data, data+n);
!       std::sort(v.begin(), v.end());
!       std::copy(v.begin(), v.end(), data);
!     }
!   }

! Fortran side:
module cpp_glue
  use iso_c_binding
  implicit none
  interface
    subroutine cpp_process(data, n) bind(C, name='cpp_process')
      import :: c_double, c_int
      real(c_double), intent(inout) :: data(*)
      integer(c_int), value :: n
    end subroutine
  end interface
end module

program use_cpp
  use cpp_glue
  implicit none
  real(c_double) :: arr(10) = [5.d0, 3.d0, 8.d0, 1.d0, 9.d0, &
                                2.d0, 7.d0, 4.d0, 6.d0, 0.d0]
  call cpp_process(arr, size(arr, kind=c_int))
  print *, arr
end program

! Build: gfortran use_cpp.f90 wrapper.cpp -lstdc++ -o app

! Fortran module (myfuncs.f90):
module myfuncs
  implicit none
contains
  subroutine compute(arr, n, result)
    integer, intent(in) :: n
    real(8), intent(in) :: arr(n)
    real(8), intent(out) :: result
    result = sum(arr**2)
  end subroutine

  function fast_sin(x) result(y)
    real(8), intent(in) :: x
    real(8) :: y
    y = sin(x)
  end function
end module

! Build Python extension:
!   f2py -c myfuncs.f90 -m myfuncs
!   (creates myfuncs.cpython-*.so)

! Python usage:
!   import numpy as np
!   import myfuncs
!   arr = np.array([1.0, 2.0, 3.0])
!   result = myfuncs.myfuncs.compute(arr)
!   y = myfuncs.myfuncs.fast_sin(1.5)

module dynamic_loader
  use iso_c_binding
  implicit none
  interface
    function dlopen(filename, flag) bind(C, name='dlopen')
      import :: c_ptr, c_char, c_int
      character(kind=c_char), dimension(*) :: filename
      integer(c_int), value :: flag
      type(c_ptr) :: dlopen
    end function

    function dlsym(handle, name) bind(C, name='dlsym')
      import :: c_ptr, c_char
      type(c_ptr), value :: handle
      character(kind=c_char), dimension(*) :: name
      type(c_ptr) :: dlsym
    end function

    function dlclose(handle) bind(C, name='dlclose')
      import :: c_ptr, c_int
      type(c_ptr), value :: handle
      integer(c_int) :: dlclose
    end function
  end interface
end module

program plugin
  use dynamic_loader
  implicit none
  type(c_ptr) :: lib, sym
  integer, parameter :: RTLD_NOW = 2

  lib = dlopen('./plugin.so' // c_null_char, RTLD_NOW)
  if (.not. c_associated(lib)) then
    print *, 'Failed to load library'
    stop 1
  end if

  sym = dlsym(lib, 'init' // c_null_char)
  ! call init via c_f_procpointer...
  print *, dlclose(lib)
end program

# CMakeLists.txt for mixed Fortran/C/C++ project:
cmake_minimum_required(VERSION 3.18)
project(mixed LANGUAGES Fortran C CXX)

set(CMAKE_Fortran_FLAGS "-O2 -fopenmp")
set(CMAKE_C_FLAGS "-O2")
set(CMAKE_CXX_FLAGS "-O2 -std=c++17")

add_library(fortran_lib STATIC mymod.f90)
add_library(c_lib STATIC helper.c)
add_library(cpp_lib STATIC wrapper.cpp)

# Mixed-language executable
add_executable(app main.f90)
target_link_libraries(app fortran_lib c_lib cpp_lib)

# Fortran needs C runtime
if(APPLE)
  target_link_libraries(app "-lc++")
else()
  target_link_libraries(app stdc++)
endif()

# Module dependency tracking (CMake handles automatically)
#   fortran_lib depends on its own .mod files
#   app depends on fortran_lib's .mod files

program safe_io
  implicit none
  integer :: u, ios
  character(256) :: msg

  ! newunit: compiler picks an unused unit number
  open(newunit=u, file='data.txt', status='old', action='read', &
       iostat=ios, iomsg=msg)
  if (ios /= 0) then
    print *, 'Open failed:', trim(msg)
    stop 1
  end if

  ! read with error handling
  read(u, *, iostat=ios, iomsg=msg) some_value
  if (ios /= 0) then
    print *, 'Read failed:', trim(msg)
  end if

  close(u)
contains
  integer function some_value()
    some_value = 0
  end function
end program

program stream_io
  implicit none
  integer :: u, i
  real, allocatable :: data(:)

  ! stream access = byte-oriented (like C fread/fwrite)
  open(newunit=u, file='data.bin', access='stream', &
       form='unformatted', status='replace')

  data = [(real(i), i=1,100)]
  write(u) data  ! binary, no record markers
  close(u)

  ! read back
  allocate(data(100))
  open(newunit=u, file='data.bin', access='stream', &
       form='unformatted', status='old', action='read')
  read(u) data
  close(u)

  ! position-based access
  open(newunit=u, file='data.bin', access='stream', &
       form='unformatted', status='old')
  read(u, pos=21) data(1)  ! read 5th real (4 bytes each)
  close(u)
end program

module person_type
  implicit none
  type :: person
    character(20) :: name
    integer :: age
    real :: height
  end type
contains
  ! custom formatted I/O via DT format
  subroutine write_person(dtv, unit, iotype, v_list, iostat, iomsg)
    class(person), intent(in) :: dtv
    integer, intent(in) :: unit
    character(*), intent(in) :: iotype
    integer, intent(in) :: v_list(:)
    integer, intent(out) :: iostat
    character(*), intent(inout) :: iomsg
    write(unit, '(a,",",i0,",",f0.2)', iostat=iostat) &
      trim(dtv%name), dtv%age, dtv%height
  end subroutine
end module

program use_dtv
  use person_type
  implicit none
  type(person) :: p = person('Alice', 30, 5.7)

  ! DT format triggers custom writer
  print "(DT)", p  ! calls write_person
end program

program namelist_demo
  implicit none
  integer :: n_iterations = 100
  real :: tolerance = 1.0e-6
  character(50) :: output_file = 'results.txt'
  logical :: verbose = .true.

  namelist /config/ n_iterations, tolerance, output_file, verbose

  ! write namelist
  open(newunit=u, file='config.nml', status='replace')
  write(u, nml=config)
  close(u)

  ! read namelist (only specified vars are updated)
  open(newunit=u, file='config.nml', status='old', action='read')
  read(u, nml=config)
  close(u)

  print *, 'Iterations:', n_iterations
end program

! config.nml format:
! &config
!   n_iterations = 500
!   tolerance = 1.0e-8
!   verbose = .false.
! /

program async_io
  implicit none
  integer :: u1, u2, ios
  real, allocatable :: a(:), b(:)

  allocate(a(1000000), b(1000000))
  a = 1.0; b = 2.0

  open(newunit=u1, file='a.bin', access='stream', form='unformatted', &
       asynchronous='yes')
  open(newunit=u2, file='b.bin', access='stream', form='unformatted', &
       asynchronous='yes')

  ! non-blocking writes
  write(u1, asynchronous='yes') a
  write(u2, asynchronous='yes') b

  ! do other work while I/O proceeds...

  ! wait for completion
  wait(u1)
  wait(u2)

  close(u1); close(u2)
end program

! gfortran debug build:
!   gfortran -O0 -g -fcheck=all -fbacktrace -ffpe-trap=invalid,zero,overflow \
!            -Wall -Wextra -Wpedantic -finit-real=nan prog.f90

! Flag explanations:
!   -O0          no optimization (easier debugging)
!   -g           debug symbols (for gdb)
!   -fcheck=all  bounds, pointer, recursion checks
!   -fbacktrace  print stack trace on error
!   -ffpe-trap   halt on FP exceptions (invalid, zero, overflow)
!   -finit-real=nan  initialize reals to NaN (catch uninitialized)
!   -finit-integer=-9999  initialize integers to sentinel
!   -Wall -Wextra  more warnings

! Runtime error message:
!   At line 42 of file prog.f90
!   Fortran runtime error: Array bound out of bounds for dimension 1

program robust_io
  implicit none
  integer :: u, ios, line_num
  character(256) :: msg, line
  real :: value

  open(newunit=u, file='data.txt', status='old', action='read', &
       iostat=ios, iomsg=msg)
  if (ios /= 0) call error_exit('Open: ' // trim(msg))

  line_num = 0
  do
    line_num = line_num + 1
    read(u, '(a)', iostat=ios) line
    if (ios < 0) exit  ! EOF
    if (ios > 0) then
      print *, 'Read error at line', line_num
      cycle
    end if

    read(line, *, iostat=ios, iomsg=msg) value
    if (ios /= 0) then
      print *, 'Parse error at line', line_num, ':', trim(msg)
      cycle
    end if
    ! process value
  end do
  close(u)
contains
  subroutine error_exit(m)
    character(*), intent(in) :: m
    print *, 'ERROR:', m
    stop 1
  end subroutine
end program

! Compile: gfortran -g -O0 prog.f90 -o prog
! Start: gdb ./prog

! Common GDB commands for Fortran:
!   (gdb) break main           break at main
!   (gdb) break prog.f90:42     break at line 42
!   (gdb) break mymod__my_sub   break at subroutine (note double underscore)
!   (gdb) run                   start program
!   (gdb) next                  step over
!   (gdb) step                  step into
!   (gdb) print arr(5)          print array element
!   (gdb) print arr(1:10)       print array slice
!   (gdb) print mat(2,3)        print 2D array element
!   (gdb) print p%name          print derived type component
!   (gdb) display x             watch variable (auto-print on stop)
!   (gdb) backtrace             show call stack
!   (gdb) info locals           show all local variables
!   (gdb) continue              resume execution

! Compile with profiling:
!   gfortran -pg -O2 prog.f90 -o prog
! Run:
!   ./prog
!   (creates gmon.out)
! Analyze:
!   gprof prog gmon.out > profile.txt
!   gprof prog gmon.out | less

! profile.txt sections:
!   - Flat profile: time per function (self + cumulative)
!   - Call graph: who called whom, how many times
!   - Index: function cross-references

! Key columns:
!   %time  percentage of total time
!   cumulative  running total
!   self    seconds in function (excluding children)
!   calls   number of calls
!   self/call  average time per call (self)

! For multi-threaded: use gprofng (modern) or perf

program benchmark
  use iso_fortran_env, only: real64, int64
  implicit none
  integer, parameter :: n = 1000000
  integer, parameter :: niter = 100
  real(real64), allocatable :: a(:), b(:), c(:)
  integer(int64) :: start, finish, rate
  real(real64) :: t_start, t_end, min_time
  integer :: i, iter

  allocate(a(n), b(n), c(n))
  a = 1.0_real64; b = 2.0_real64

  ! warmup (cache, JIT-like effects)
  do i = 1, n
    c(i) = a(i) + b(i)
  end do

  min_time = huge(min_time)
  call system_clock(count_rate=rate)
  do iter = 1, niter
    call system_clock(start)
    do i = 1, n
      c(i) = a(i) + b(i)
    end do
    call system_clock(finish)
    t_start = real(start, real64) / rate
    t_end = real(finish, real64) / rate
    min_time = min(min_time, t_end - t_start)
  end do

  print *, 'Best time:', min_time * 1e6, 'us'
  print *, 'Throughput:', real(n) / min_time / 1e9, 'G elem/s'
end program

program modern
    implicit none
    integer :: i
    do i = 1, 10
        print *, "Value:", i
    end do
end program modern
! Free form: no column restrictions
! Comments start with !

module math_utils
    implicit none
    contains
    function square(x) result(y)
        real, intent(in) :: x
        real :: y
        y = x * x
    end function square
end module math_utils

program test
    use math_utils
    print *, square(3.0)  ! 9.0
end program test

type :: Point
    real :: x, y
end type Point
type(Point) :: p
p = Point(1.0, 2.0)  ! Constructor
p%x = 3.0  ! Component access
print *, p%x, p%y
! Type-bound procedures (OOP)
type :: Circle
    real :: radius
contains
    procedure :: area => circle_area
end type

subroutine process(input, output, inout)
    integer, intent(in) :: input    ! Read-only
    integer, intent(out) :: output  ! Write-only
    integer, intent(inout) :: inout ! Read-write
    output = input * 2
    inout = inout + 1
end subroutine

pure function square(x) result(y)
    real, intent(in) :: x
    real :: y
    y = x * x
end function
! Elemental: works on scalars and arrays
elemental function double_it(x) result(y)
    real, intent(in) :: x
    real :: y
    y = 2.0 * x
end function
! double_it([1,2,3]) returns [2,4,6]

! Static arrays
real :: a(10)           ! 1D, indices 1-10
real :: b(0:9)          ! 1D, indices 0-9
real :: c(3, 4)         ! 2D, 3 rows x 4 cols
! Allocatable (dynamic)
real, allocatable :: d(:)
allocate(d(100))        ! Allocate
deallocate(d)           ! Free
! Array constants
integer :: nums(5) = [1, 2, 3, 4, 5]

real :: a(5) = [1, 2, 3, 4, 5]
real :: b(5)
b = a * 2          ! Element-wise: [2,4,6,8,10]
b = a + 1          ! [2,3,4,5,6]
print *, sum(a)    ! 15
print *, maxval(a) ! 5
print *, any(a > 3)  ! .true.
print *, count(a > 2) ! 3

real :: a(10) = [(i, i=1,10)]
print *, a(3:7)     ! Elements 3 to 7
print *, a(::2)     ! Every other: [1,3,5,7,9]
print *, a(2:8:2)   ! Stride 2: [2,4,6,8]
real :: m(3,3)
m(:, 2)             ! Second column
m(1, :)             ! First row

real :: a(5) = [1, -2, 3, -4, 5]
where (a > 0)
    a = a * 2       ! Double positives
elsewhere
    a = 0           ! Zero negatives
end where
! Result: [2, 0, 6, 0, 10]
! Equivalent to a loop with if

real, allocatable :: matrix(:,:)
integer :: n
n = 100
allocate(matrix(n, n))
matrix = 0.0  ! Initialize all to 0
! ... use matrix ...
deallocate(matrix)
! Automatic deallocation at end of scope
! (Fortran 2003+)

integer :: i = 42
real :: x = 3.14159
print "(I5, F8.3)", i, x    ! "   42    3.142"
print "(A, I0)", "Count=", i  ! "Count=42"
write(*, "(3F6.2)") 1.0, 2.0, 3.0
! Format specifiers:
! I: integer, F: float, E: exponential
! A: string, X: space, /: newline

integer :: unit, ios
open(newunit=unit, file="data.txt", status="old", action="read")
do
    read(unit, *, iostat=ios) value
    if (ios /= 0) exit
    print *, value
end do
close(unit)
! status: old, new, replace, scratch
! action: read, write, readwrite

integer :: n = 10
real :: x = 3.14
namelist /config/ n, x
! Write namelist
open(1, file="config.nml")
write(1, nml=config)
close(1)
! Read namelist
open(1, file="config.nml")
read(1, nml=config)
close(1)

character(20) :: str
integer :: num = 42
! Integer to string
write(str, "(I0)") num
print *, trim(str)  ! "42"
! String to integer
read(str, *) num
! Internal files use character variables as units

! Unformatted (binary) I/O
open(1, file="data.bin", form="unformatted", &
     access="stream")
write(1) array  ! No format, raw bytes
read(1) array2
close(1)
! Faster than formatted I/O
! Smaller file size
! Not portable across architectures

!$omp parallel do private(i) reduction(+:sum)
do i = 1, n
    sum = sum + a(i) * b(i)
end do
!$omp end parallel do
! Compile: gfortran -fopenmp program.f90
! Environment: OMP_NUM_THREADS=4

program coarray_example
    implicit none
    integer :: me[*]  ! Coarray: one per image
    me = this_image()
    sync all  ! Barrier
    if (this_image() == 1) then
        print *, "Image 2 has:", me[2]  ! Remote access
    end if
end program
! Compile: gfortran -fcoarray=single program.f90

program mpi_example
    use mpi
    integer :: rank, size, ierr
    call mpi_init(ierr)
    call mpi_comm_rank(MPI_COMM_WORLD, rank, ierr)
    call mpi_comm_size(MPI_COMM_WORLD, size, ierr)
    print *, "I am rank", rank, "of", size
    call mpi_finalize(ierr)
end program
! Compile: mpifort program.f90

do concurrent (i = 1:n) local(tmp)
    tmp = expensive_compute(a(i))
    b(i) = tmp * 2
end do
! Tells compiler iterations are independent
! Can be parallelized automatically
! local: private variable per iteration

!$omp parallel do reduction(+:total)
do i = 1, n
    total = total + a(i)
end do
!$omp end parallel do
! Common reductions: +, *, max, min, .and., .or.
! Each thread has a private copy
! Combined at the end

! Matrix multiplication
do i = 1, n
    do j = 1, n
        c(i,j) = sum(a(i,:) * b(:,j))
    end do
end do
! Or use matmul intrinsic
c = matmul(a, b)
! Dot product
dot = dot_product(a, b)
! Transpose
at = transpose(a)

! Euler method: dy/dt = f(t, y)
subroutine euler(f, t0, y0, dt, n, t, y)
    interface
        real function f(t, y)
            real, intent(in) :: t, y
        end function
    end interface
    real, intent(in) :: t0, y0, dt
    integer, intent(in) :: n
    real, intent(out) :: t(n+1), y(n+1)
    integer :: i
    t(1) = t0; y(1) = y0
    do i = 1, n
        t(i+1) = t(i) + dt
        y(i+1) = y(i) + dt * f(t(i), y(i))
    end do
end subroutine

call random_seed()  ! Seed from system
call random_number(x)  ! x in [0, 1)
! Array of randoms
real :: arr(100)
call random_number(arr)
! Integer in range [1, 6]
integer :: dice
call random_number(r)
dice = int(r * 6) + 1
! Normal distribution (Box-Muller)
call random_number(u1)
call random_number(u2)
z = sqrt(-2*log(u1)) * cos(2*PI*u2)

function interp(x, xs, ys) result(y)
    real, intent(in) :: x, xs(:), ys(:)
    real :: y
    integer :: i
    ! Find interval
    i = 1
    do while (i < size(xs) .and. x > xs(i+1))
        i = i + 1
    end do
    ! Linear interpolation
    y = ys(i) + (ys(i+1) - ys(i)) * &
        (x - xs(i)) / (xs(i+1) - xs(i))
end function

! Trapezoidal rule
function trapezoid(f, a, b, n) result(integral)
    interface
        real function f(x)
            real, intent(in) :: x
        end function
    end interface
    real, intent(in) :: a, b
    integer, intent(in) :: n
    real :: integral, h
    integer :: i
    h = (b - a) / n
    integral = (f(a) + f(b)) / 2
    do i = 1, n-1
        integral = integral + f(a + i*h)
    end do
    integral = integral * h
end function

! Fortran arrays are 1-based by default
real :: a(10)  ! Indices 1 to 10
! a(0) = 1.0  ! Error: out of bounds
a(1) = 1.0    ! OK
! Custom bounds
real :: b(0:9)  ! Indices 0 to 9
b(0) = 1.0    ! OK

! BAD: implicit typing (Fortran 77 style)
program bad
    ! i-n start with integer by default
    i = 1      ! integer (implicit)
    x = 3.14   ! real (implicit)
end program
! GOOD: explicit typing
program good
    implicit none  ! Force explicit declaration
    integer :: i
    real :: x
end program

! Fortran is column-major (first index varies fastest)
real :: m(3, 3)
! Efficient: iterate over first index in inner loop
do j = 1, 3
    do i = 1, 3
        m(i, j) = 0.0  ! Cache-friendly
    end do
end do
! Inefficient: row-major access
do i = 1, 3
    do j = 1, 3
        m(i, j) = 0.0  ! Cache misses
    end do
end do

subroutine modify(x)
    integer, intent(inout) :: x
    x = 99  ! Modifies the caller variable
end subroutine
! Fortran passes by reference by default
! intent(in) prevents modification
! Without intent, modification is allowed (dangerous)

! Single precision (default)
real :: x = 3.14159  ! ~7 digits
! Double precision
real(kind=8) :: y = 3.14159d0  ! ~15 digits
! Or use kind parameter
integer, parameter :: dp = kind(1.0d0)
real(dp) :: z = 3.14159_dp
! Always use _dp or d0 for double literals