Expressions and Assignment

module heatflow_kinds
  use, intrinsic :: iso_fortran_env, only: int32, int64, real32, real64
  implicit none
  private
  public :: sp, dp, ip, lp
  integer, parameter :: sp = real32
  integer, parameter :: dp = real64
  integer, parameter :: ip = int32
  integer, parameter :: lp = int64
end module heatflow_kinds

program precedence_demo
  use heatflow_kinds, only: dp, ip
  implicit none
  integer(ip) :: two, three
  real(dp) :: a
  logical :: p, q, r
  character(len=:), allocatable :: label

  two = 2_ip
  three = 3_ip
  a = 3.0_dp
  p = .true.
  q = .false.
  r = .false.
  label = 'ab'

  ! Numeric: ** binds tighter than unary minus and is right-associative.
  print '(a, i0)', '-two**2        = ', -two**2
  print '(a, i0)', '(-two)**2      = ', (-two)**2
  print '(a, i0)', 'two**three**2  = ', two**three**2
  print '(a, i0)', '(two**three)**2 = ', (two**three)**2
  print '(a, f0.1)', '-a*a           = ', -a*a
  ! Numeric operators bind tighter than relational ones.
  print '(a, l1)', 'two + 1 == three = ', two + 1_ip == three
  ! Concatenation binds tighter than relational ones.
  print '(a, l1)', 'label // "c" == "abc" = ', label // 'c' == 'abc'
  ! .not. is below relational operators, so it applies to the whole comparison.
  print '(a, l1)', '.not. two == three = ', .not. two == three
  ! .and. binds tighter than .or., which binds tighter than .eqv. / .neqv.
  print '(a, l1)', 'p .or. q .and. r   = ', p .or. q .and. r
  print '(a, l1)', '(p .or. q) .and. r = ', (p .or. q) .and. r
  print '(a, l1)', 'q .eqv. r .or. p   = ', q .eqv. r .or. p
  print '(a, l1)', '(q .eqv. r) .or. p = ', (q .eqv. r) .or. p
end program precedence_demo

-two**2        = -4
(-two)**2      = 4
two**three**2  = 512
(two**three)**2 = 64
-a*a           = -9.0
two + 1 == three = T
label // "c" == "abc" = T
.not. two == three = T
p .or. q .and. r   = T
(p .or. q) .and. r = F
q .eqv. r .or. p   = F
(q .eqv. r) .or. p = T

! REJECTED on purpose: neither statement is valid Fortran.
program rejected_operator_usage
  implicit none
  logical :: a, b
  real :: x, y
  a = .true.
  b = .false.
  x = 2.0
  y = 3.0
  print *, a == b
  print *, x ** -y
end program rejected_operator_usage

16_rejected_operator_usage.f90:11:18:

   11 |   print *, x ** -y
      |                  1
Error: Extension: Unary operator following arithmetic operator (use parentheses) at (1)
16_rejected_operator_usage.f90:10:10:

   10 |   print *, a == b
      |          1
Error: Logicals at (1) must be compared with .eqv. instead of ==

error: Semantic errors in 16_rejected_operator_usage.f90
16_rejected_operator_usage.f90:10:12: error: LOGICAL operands must be compared using .EQV. or .NEQV.
    print *, a == b
             ^^^^^^

module heatflow_kinds
  use, intrinsic :: iso_fortran_env, only: int32, int64, real32, real64
  implicit none
  private
  public :: sp, dp, ip, lp
  integer, parameter :: sp = real32
  integer, parameter :: dp = real64
  integer, parameter :: ip = int32
  integer, parameter :: lp = int64
end module heatflow_kinds

program numeric_operations
  use heatflow_kinds, only: sp, dp, ip
  implicit none
  integer(ip) :: nx, ny, i, j
  real(dp) :: lx, dx, x
  real(sp) :: tenth_sp
  real(dp) :: tenth_dp

  ! Integer division truncates toward zero (10.1.5.2.2); mod follows the
  ! dividend, modulo follows the divisor.
  i = -8_ip
  j = 3_ip
  print '(a, i0)', 'i / j          = ', i/j
  print '(a, i0)', 'mod(i, j)      = ', mod(i, j)
  print '(a, i0)', 'modulo(i, j)   = ', modulo(i, j)
  print '(a, i0)', 'j ** (-2)      = ', j**(-2_ip)

  ! Mixed mode: the integer operand is converted to the real kind of the
  ! other operand. Convert before dividing to avoid integer division.
  nx = 5_ip
  ny = 2_ip
  lx = 1.0_dp
  dx = lx/real(nx - 1_ip, dp)
  print '(a, i0)', 'nx / ny        = ', nx/ny
  print '(a, f6.4)', 'real(nx, dp)/ny = ', real(nx, dp)/ny
  print '(a, f6.4)', 'dx = lx/(nx-1)  = ', dx

  ! Operators of equal precedence combine left to right, so the integer
  ! division in the first form happens before the multiplication.
  x = 2.5_dp
  print '(a, f6.4)', 'nx / ny * x    = ', nx/ny*x
  print '(a, f6.4)', 'nx * x / ny    = ', nx*x/ny

  ! A default real constant is single precision; its error survives
  ! assignment to a double precision variable.
  tenth_sp = 0.1
  tenth_dp = real(tenth_sp, dp)
  print '(a, es24.16)', 'real(0.1, dp)  = ', tenth_dp
  tenth_dp = 0.1_dp
  print '(a, es24.16)', '0.1_dp         = ', tenth_dp
end program numeric_operations

i / j          = -2
mod(i, j)      = -2
modulo(i, j)   = 1
j ** (-2)      = 0
nx / ny        = 2
real(nx, dp)/ny = 2.5000
dx = lx/(nx-1)  = 0.2500
nx / ny * x    = 5.0000
nx * x / ny    = 6.2500
real(0.1, dp)  =   1.0000000149011612E-01
0.1_dp         =   1.0000000000000001E-01

! NON-CONFORMING on purpose: division by zero is not defined by the
! arithmetic of the processor (10.1.5.2.4), and -fcheck=all does not catch it.
program integer_division_by_zero
  implicit none
  integer :: i, j
  i = 1
  j = 0
  print '(a, i0)', 'i / j = ', i/j
end program integer_division_by_zero

i / j = 0

module heatflow_kinds
  use, intrinsic :: iso_fortran_env, only: int32, int64, real32, real64
  implicit none
  private
  public :: sp, dp, ip, lp
  integer, parameter :: sp = real32
  integer, parameter :: dp = real64
  integer, parameter :: ip = int32
  integer, parameter :: lp = int64
end module heatflow_kinds

program relational_and_logical
  use, intrinsic :: ieee_arithmetic, only: ieee_value, ieee_quiet_nan
  use heatflow_kinds, only: dp
  implicit none
  real(dp) :: x, y, nan, residual, tolerance
  logical :: converged, hot, stable
  character(len=8) :: solver

  ! Character comparison pads the shorter operand with blanks (10.1.5.5.1).
  solver = 'jacobi'
  print '(a, l1)', "solver == 'jacobi'          = ", solver == 'jacobi'
  print '(a, l1)', "'abc' == 'abc   '           = ", 'abc' == 'abc   '
  print '(a, l1)', "'Jacobi' < 'jacobi'         = ", 'Jacobi' < 'jacobi'
  print '(a, l1)', "lgt('b', 'a')               = ", lgt('b', 'a')

  ! Logical values are compared with .eqv. and .neqv., never with ==.
  hot = .true.
  stable = .false.
  print '(a, l1)', 'hot .eqv. stable            = ', hot .eqv. stable
  print '(a, l1)', 'hot .neqv. stable           = ', hot .neqv. stable

  ! Real values: compare a difference with a tolerance instead of testing
  ! equality (gfortran -Wextra warns about == on reals).
  x = 0.1_dp + 0.2_dp
  y = 0.3_dp
  print '(a, es10.3)', 'abs(x - y)                  = ', abs(x - y)
  tolerance = 1.0e-12_dp
  residual = abs(x - y)
  converged = residual <= tolerance
  print '(a, l1)', 'residual <= tolerance       = ', converged

  ! With a NaN operand, .not. (a < b) and a >= b are different expressions.
  nan = ieee_value(nan, ieee_quiet_nan)
  x = 1.0_dp
  print '(a, l1)', '.not. (nan < x)             = ', .not. (nan < x)
  print '(a, l1)', 'nan >= x                    = ', nan >= x
end program relational_and_logical

solver == 'jacobi'          = T
'abc' == 'abc   '           = T
'Jacobi' < 'jacobi'         = T
lgt('b', 'a')               = T
hot .eqv. stable            = F
hot .neqv. stable           = T
abs(x - y)                  =  5.551E-17
residual <= tolerance       = T
.not. (nan < x)             = T
nan >= x                    = F

module heatflow_kinds
  use, intrinsic :: iso_fortran_env, only: int32, int64, real32, real64
  implicit none
  private
  public :: sp, dp, ip, lp
  integer, parameter :: sp = real32
  integer, parameter :: dp = real64
  integer, parameter :: ip = int32
  integer, parameter :: lp = int64
end module heatflow_kinds

module temperature_demo
  use heatflow_kinds, only: dp
  implicit none
  private
  public :: temperature_t, operator(+), operator(.celsius.), operator(.diff.), assignment(=)

  ! A temperature stored in kelvin; the type keeps scales from being mixed.
  type :: temperature_t
    real(dp) :: kelvin = 0.0_dp
  end type temperature_t

  interface operator(+)
    module procedure add_delta
  end interface

  interface operator(.celsius.)
    module procedure to_celsius
  end interface

  interface operator(.diff.)
    module procedure difference
  end interface

  interface assignment(=)
    module procedure set_from_celsius
  end interface

contains

  pure function add_delta(t, delta) result(r)
    type(temperature_t), intent(in) :: t
    real(dp), intent(in) :: delta
    type(temperature_t) :: r
    r%kelvin = t%kelvin + delta
  end function add_delta

  pure function to_celsius(t) result(c)
    type(temperature_t), intent(in) :: t
    real(dp) :: c
    c = t%kelvin - 273.15_dp
  end function to_celsius

  pure function difference(a, b) result(d)
    type(temperature_t), intent(in) :: a, b
    real(dp) :: d
    d = a%kelvin - b%kelvin
  end function difference

  ! t = 21.5_dp stores a Celsius value as kelvin.
  pure subroutine set_from_celsius(t, celsius)
    type(temperature_t), intent(out) :: t
    real(dp), intent(in) :: celsius
    t%kelvin = celsius + 273.15_dp
  end subroutine set_from_celsius

end module temperature_demo

program defined_operators_and_assignment
  use heatflow_kinds, only: dp
  use temperature_demo, only: temperature_t, operator(+), operator(.celsius.), &
    operator(.diff.), assignment(=)
  implicit none
  type(temperature_t) :: ambient, hot, warmer

  ambient = 20.0_dp
  hot = 80.0_dp
  warmer = ambient + 5.0_dp
  print '(a, f7.2)', 'ambient in Celsius:      ', .celsius. ambient
  print '(a, f7.2)', 'warmer in Celsius:       ', .celsius. warmer
  ! A defined binary operator has the lowest precedence, so both sums are
  ! formed first: (hot + 0) .diff. (ambient + 0).
  print '(a, f7.2)', 'hot .diff. ambient:      ', hot + 0.0_dp .diff. ambient + 0.0_dp
  ! A defined unary operator has the highest precedence: it applies to
  ! "ambient", not to the sum.
  print '(a, f7.2)', '.celsius. ambient + 1.0: ', .celsius. ambient + 1.0_dp
end program defined_operators_and_assignment

ambient in Celsius:        20.00
warmer in Celsius:         25.00
hot .diff. ambient:        60.00
.celsius. ambient + 1.0:   21.00

( scalar-logical-expr ? expr [ : scalar-logical-expr ? expr ]... : expr )

module heatflow_kinds
  use, intrinsic :: iso_fortran_env, only: int32, int64, real32, real64
  implicit none
  private
  public :: sp, dp, ip, lp
  integer, parameter :: sp = real32
  integer, parameter :: dp = real64
  integer, parameter :: ip = int32
  integer, parameter :: lp = int64
end module heatflow_kinds

module guard_demo
  use heatflow_kinds, only: dp, ip
  implicit none
  private
  public :: positive_at, sample, describe

contains

  ! Safe guard: the subscript is only evaluated when the condition is true.
  pure function positive_at(a, i) result(ok)
    real(dp), intent(in) :: a(:)
    integer(ip), intent(in) :: i
    logical :: ok
    ok = (i >= 1_ip .and. i <= size(a, kind=ip) ? a(i) > 0.0_dp : .false.)
  end function positive_at

  ! Optional argument with a default, without a second variable.
  pure function sample(a, i, fill) result(v)
    real(dp), intent(in) :: a(:)
    integer(ip), intent(in) :: i
    real(dp), intent(in), optional :: fill
    real(dp) :: v
    v = (i >= 1_ip .and. i <= size(a, kind=ip) ? a(i) : present(fill) ? fill : 0.0_dp)
  end function sample

  ! The branches may have different character lengths; the result takes the
  ! length of the chosen one.
  function describe(residual, tolerance) result(text)
    real(dp), intent(in) :: residual, tolerance
    character(len=:), allocatable :: text
    text = (abs(residual) <= tolerance ? 'ok' : 'did not converge')
  end function describe

end module guard_demo

program conditional_expressions
  use heatflow_kinds, only: dp, ip
  use guard_demo, only: positive_at, sample, describe
  implicit none
  real(dp) :: a(3)
  integer(ip) :: i

  a = [1.0_dp, -2.0_dp, 3.0_dp]
  do i = 0_ip, 4_ip
    print '(a, i0, a, l1, a, f6.2, a, f6.2)', 'i = ', i, '  positive: ', positive_at(a, i), &
      '  sample: ', sample(a, i), '  sample(fill=-1): ', sample(a, i, -1.0_dp)
  end do

  print '(a)', describe(2.0e-5_dp, 1.0e-4_dp)
  print '(a)', describe(3.0e-2_dp, 1.0e-4_dp)

  ! Chain of conditions, evaluated in order.
  print '(a)', (a(1) > 5.0_dp ? 'big' : a(1) > 0.5_dp ? 'medium' : 'small')
end program conditional_expressions

i = 0  positive: F  sample:   0.00  sample(fill=-1):  -1.00
i = 1  positive: T  sample:   1.00  sample(fill=-1):   1.00
i = 2  positive: F  sample:  -2.00  sample(fill=-1):  -2.00
i = 3  positive: T  sample:   3.00  sample(fill=-1):   3.00
i = 4  positive: F  sample:   0.00  sample(fill=-1):  -1.00
ok
did not converge
medium

! REJECTED on purpose: the branches of a conditional expression must have
! the same declared type and kind (C1004).
program conditional_mixed_types
  implicit none
  integer :: i = 1
  print *, (i > 0 ? 1 : 2.0)
end program conditional_mixed_types

17_conditional_expression_rejected.f90:6:20-23:

    6 |   print *, (i > 0 ? 1 : 2.0)
      |                    1  2
Error: expr at (1) and expr at (2) in conditional expression must have the same declared type

error: Semantic errors in 17_conditional_expression_rejected.f90
17_conditional_expression_rejected.f90:6:12: error: All values in conditional expression must have the same type and kind; have INTEGER(4) and REAL(4)
    print *, (i > 0 ? 1 : 2.0)
             ^^^^^^^^^^^^^^^^^

module heatflow_kinds
  use, intrinsic :: iso_fortran_env, only: int32, int64, real32, real64
  implicit none
  private
  public :: sp, dp, ip, lp
  integer, parameter :: sp = real32
  integer, parameter :: dp = real64
  integer, parameter :: ip = int32
  integer, parameter :: lp = int64
end module heatflow_kinds

module field_tools
  use heatflow_kinds, only: dp, ip
  implicit none
  private
  public :: coarsen, report

contains

  ! Keep every second point in each direction.
  pure function coarsen(u) result(c)
    real(dp), intent(in) :: u(:, :)
    real(dp), allocatable :: c(:, :)
    c = u(1:size(u, 1):2, 1:size(u, 2):2)
  end function coarsen

  subroutine report(name, u)
    character(len=*), intent(in) :: name
    real(dp), allocatable, intent(in) :: u(:, :)
    if (allocated(u)) then
      print '(a, a, i0, a, i0, a, i0, a, i0)', name, ': shape ', size(u, 1), ' x ', size(u, 2), &
        ', lbound ', lbound(u, 1), ', ', lbound(u, 2)
    else
      print '(a, a)', name, ': not allocated'
    end if
  end subroutine report

end module field_tools

program allocatable_assignment
  use heatflow_kinds, only: dp, ip
  use field_tools, only: coarsen, report
  implicit none
  real(dp), allocatable :: u(:, :), unew(:, :)
  real(dp) :: fine(9, 9)
  real(dp) :: halo(0:4, 0:3)
  integer(ip) :: i

  fine = reshape([(real(i, dp), i=1_ip, 81_ip)], shape(fine))

  call report('u', u)
  ! Unallocated target: allocated with the shape of the right-hand side.
  u = fine
  call report('u = fine', u)
  ! Same shape: no reallocation, just the copy.
  unew = u
  call report('unew = u', unew)
  ! Different shape: u is deallocated and allocated again.
  u = coarsen(fine)
  call report('u = coarsen(fine)', u)
  ! A scalar is broadcast to an allocated array; the shape is kept.
  u = 0.0_dp
  call report('u = 0', u)
  ! Section on the left: no reallocation, the shapes must conform.
  u(:, :) = coarsen(fine)
  call report('u(:,:) = ...', u)
  ! A section or a function result has lower bounds 1 ...
  u = fine(3:5, 4:6)
  call report('u = fine(3:5, 4:6)', u)
  ! ... but a whole array keeps its declared lower bounds, and the
  ! reallocated target takes them over.
  halo = 1.0_dp
  u = halo
  call report('u = halo', u)
  deallocate (u)
  call report('deallocate(u)', u)
  print '(a, f6.1)', 'unew(9, 9) = ', unew(9, 9)
end program allocatable_assignment

u: not allocated
u = fine: shape 9 x 9, lbound 1, 1
unew = u: shape 9 x 9, lbound 1, 1
u = coarsen(fine): shape 5 x 5, lbound 1, 1
u = 0: shape 5 x 5, lbound 1, 1
u(:,:) = ...: shape 5 x 5, lbound 1, 1
u = fine(3:5, 4:6): shape 3 x 3, lbound 1, 1
u = halo: shape 5 x 4, lbound 0, 0
deallocate(u): not allocated
unew(9, 9) =   81.0

! NON-CONFORMING on purpose: a scalar cannot define the shape of an
! unallocated array target (10.2.1.3).
program scalar_into_unallocated
  implicit none
  real, allocatable :: u(:, :)
  u = 0.0
  print '(a, l1)', 'allocated(u) = ', allocated(u)
end program scalar_into_unallocated

08a_scalar_into_unallocated_nonconforming.f90:6:9:
Warning: 'u.dim[0].ubound' may be used uninitialized [-Wmaybe-uninitialized]
(three more warnings of the same kind)
allocated(u) = F

fatal Fortran runtime error(08a_scalar_into_unallocated_nonconforming.f90:6): Assign: mismatched ranks (2 != 0) in assignment to unallocated allocatable
Aborted

! NON-CONFORMING on purpose: u(:, :) is a section, so no reallocation
! happens and the shapes must conform (10.2.1.2 item 3).
program section_shape_mismatch
  implicit none
  real, allocatable :: u(:, :), v(:, :)
  allocate (u(4, 4), v(3, 3))
  u = 0.0
  v = 1.0
  u(:, :) = v
  print '(a, f4.1)', 'u(1, 1) = ', u(1, 1)
end program section_shape_mismatch

At line 9 of file 08b_section_shape_mismatch_nonconforming.f90
Fortran runtime error: Array bound mismatch for dimension 1 of array 'u' (4/3)

fatal Fortran runtime error(08b_section_shape_mismatch_nonconforming.f90:9): Assign: mismatching element counts in array assignment (to 16, from 9)
Aborted

module heatflow_kinds
  use, intrinsic :: iso_fortran_env, only: int32, int64, real32, real64
  implicit none
  private
  public :: sp, dp, ip, lp
  integer, parameter :: sp = real32
  integer, parameter :: dp = real64
  integer, parameter :: ip = int32
  integer, parameter :: lp = int64
end module heatflow_kinds

program array_assignment_and_where
  use heatflow_kinds, only: dp, ip
  implicit none
  integer(ip), parameter :: nx = 5_ip, ny = 5_ip
  real(dp), parameter :: dx = 0.25_dp, dy = 0.25_dp, alpha = 1.0_dp, dt = 0.01_dp
  real(dp) :: u(nx, ny), unew(nx, ny)
  real(dp) :: line(6)
  integer(ip) :: j

  ! A scalar is broadcast to every element.
  u = 0.0_dp
  u(3, 3) = 1.0_dp

  ! Whole-array arithmetic on conformable sections: one FTCS step of the heat
  ! equation. Every right-hand side is evaluated before unew is changed.
  unew = u
  unew(2:nx - 1, 2:ny - 1) = u(2:nx - 1, 2:ny - 1) + alpha*dt*( &
    (u(3:nx, 2:ny - 1) - 2.0_dp*u(2:nx - 1, 2:ny - 1) + u(1:nx - 2, 2:ny - 1))/dx**2 + &
    (u(2:nx - 1, 3:ny) - 2.0_dp*u(2:nx - 1, 2:ny - 1) + u(2:nx - 1, 1:ny - 2))/dy**2)
  do j = 1_ip, ny
    print '(5f8.4)', unew(:, j)
  end do

  ! Overlapping sections: the whole right-hand side is evaluated first.
  line = [1.0_dp, 2.0_dp, 3.0_dp, 4.0_dp, 5.0_dp, 6.0_dp]
  line(1:6) = line(6:1:-1)
  print '(6f5.1)', line

  ! Masked assignment: clip temperatures, element by element.
  where (unew > 0.3_dp)
    unew = 0.3_dp
  elsewhere (unew < 0.01_dp)
    unew = 0.0_dp
  end where
  print '(a, f6.3, a, f6.3)', 'max = ', maxval(unew), '  min = ', minval(unew)
  print '(a, i0)', 'cells above zero = ', count(unew > 0.0_dp)
end program array_assignment_and_where

  0.0000  0.0000  0.0000  0.0000  0.0000
  0.0000  0.0000  0.1600  0.0000  0.0000
  0.0000  0.1600  0.3600  0.1600  0.0000
  0.0000  0.0000  0.1600  0.0000  0.0000
  0.0000  0.0000  0.0000  0.0000  0.0000
  6.0  5.0  4.0  3.0  2.0  1.0
max =  0.300  min =  0.000
cells above zero = 5

! REJECTED on purpose: the shapes of a and b are not conformable.
program nc
implicit none
real :: a(4), b(3)
a = 0.0
b = 1.0
a = b
print *, a
end program

10_nonconformable_arrays_rejected.f90:7:1:
    7 | a = b
      | 1~
Error: Different shape for array assignment at (1) on dimension 1 (4 and 3)

module heatflow_kinds
  use, intrinsic :: iso_fortran_env, only: int32, int64, real32, real64
  implicit none
  private
  public :: sp, dp, ip, lp
  integer, parameter :: sp = real32
  integer, parameter :: dp = real64
  integer, parameter :: ip = int32
  integer, parameter :: lp = int64
end module heatflow_kinds

program derived_type_assignment
  use heatflow_kinds, only: dp, ip
  implicit none

  type :: probe_t
    integer(ip) :: id = 0_ip
    real(dp), allocatable :: readings(:)    ! owned by the object
    real(dp), pointer :: shared(:) => null() ! only associated, never owned
  end type probe_t

  real(dp), target :: log_buffer(3)
  type(probe_t) :: a, b

  log_buffer = [10.0_dp, 20.0_dp, 30.0_dp]
  a%id = 7_ip
  a%readings = [0.1_dp, 0.2_dp]
  a%shared => log_buffer

  ! Intrinsic assignment: ordinary components are copied, allocatable
  ! components are deep-copied (and reallocated if needed), pointer
  ! components are pointer-assigned (10.2.1.3).
  b = a

  a%readings(1) = -1.0_dp
  a%shared(1) = -10.0_dp

  print '(a, i0)', 'b%id            = ', b%id
  print '(a, 2f6.1)', 'b%readings      =', b%readings
  print '(a, 3f6.1)', 'b%shared        =', b%shared
  print '(a, l1)', 'same pointer target: ', associated(b%shared, a%shared)
  b%readings(2) = 5.0_dp
  print '(a, 2f6.1)', 'a%readings      =', a%readings
end program derived_type_assignment

b%id            = 7
b%readings      =   0.1   0.2
b%shared        = -10.0  20.0  30.0
same pointer target: T
a%readings      =  -1.0   0.2

module heatflow_kinds
  use, intrinsic :: iso_fortran_env, only: int32, int64, real32, real64
  implicit none
  private
  public :: sp, dp, ip, lp
  integer, parameter :: sp = real32
  integer, parameter :: dp = real64
  integer, parameter :: ip = int32
  integer, parameter :: lp = int64
end module heatflow_kinds

module kernels_demo
  use heatflow_kinds, only: dp
  implicit none
  private
  public :: scale_by_two, cool

contains

  pure function scale_by_two(t) result(r)
    real(dp), intent(in) :: t
    real(dp) :: r
    r = 2.0_dp*t
  end function scale_by_two

  pure function cool(t) result(r)
    real(dp), intent(in) :: t
    real(dp) :: r
    r = 0.5_dp*t
  end function cool

end module kernels_demo

program pointer_assignment
  use heatflow_kinds, only: dp, ip
  use kernels_demo, only: scale_by_two, cool
  implicit none
  real(dp), target :: u(4, 3)
  real(dp), target :: flat(12)
  real(dp), pointer :: row(:), interior(:, :), matrix(:, :), p
  procedure(scale_by_two), pointer :: update
  integer(ip) :: i

  u = reshape([(real(i, dp), i=1_ip, 12_ip)], shape(u))

  ! Pointer assignment (=>) associates; it does not copy.
  row => u(2, :)
  row(1) = -1.0_dp
  print '(a, f5.1)', 'u(2, 1) after row(1) = -1: ', u(2, 1)

  ! Ordinary assignment through a pointer assigns to its target.
  p => u(4, 3)
  p = 99.0_dp
  print '(a, f5.1)', 'u(4, 3) after p = 99:      ', u(4, 3)

  ! A pointer to a section has lower bounds 1 unless you give them.
  interior => u(2:3, 2:3)
  print '(a, 2i3, a, 2i3)', 'lbound(interior) =', lbound(interior), '   shape =', shape(interior)
  interior(0:, 0:) => u(2:3, 2:3)
  print '(a, 2i3)', 'lbound with (0:, 0:) =', lbound(interior)

  ! Bounds remapping: a rank-2 view of rank-1 storage.
  flat = [(real(i, dp), i=1_ip, 12_ip)]
  matrix(1:4, 1:3) => flat
  print '(a, f5.1)', 'matrix(2, 3) = ', matrix(2, 3)

  ! Pointer status.
  print '(a, l1)', 'associated(row, u(2, :)) = ', associated(row, u(2, :))
  nullify (row)
  print '(a, l1)', 'associated(row) after nullify = ', associated(row)

  ! Procedure pointer assignment.
  update => scale_by_two
  print '(a, f5.1)', 'update(3) = ', update(3.0_dp)
  update => cool
  print '(a, f5.1)', 'update(3) = ', update(3.0_dp)
end program pointer_assignment

u(2, 1) after row(1) = -1:  -1.0
u(4, 3) after p = 99:       99.0
lbound(interior) =  1  1   shape =  2  2
lbound with (0:, 0:) =  0  0
matrix(2, 3) =  10.0
associated(row, u(2, :)) = T
associated(row) after nullify = F
update(3) =   6.0
update(3) =   1.5