Replacement for obsolete "statement functions" in Fortran 95 and later?

Question

Original Question

As I understand, as of Fortran 95 Statement functions have been declared obsolete in favor of internal functions. However, internal functions do not cover all use-cases, in particular when using statement functions as abbreviations with implied arguments for making the implementation of long formulas more readable. Is there any adequate replacement for this usecase?

Example

E.g., lets say we have an relation

         Cᵢⱼ = ∫dx₁∫dx₂∫dx₃∫dy₁∫dy₂∫dy₃ (AᵢⱼBᵢⱼ +BᵢⱼAⱼᵢ)

E.g. compare an implementation using directly the internal representation of A,B with the index-order mandated by prior code, resulting in

do a1=1,NGRID1; do a2=1,NGRID2; ...; do aN=1,NGRIDN
  x(i,j) = x(i,j) &
      & + ARR_A(x1,x2,x3,i,y1,y2,y3,j)*ARR_B(x1,x2,x3,i,y1,y2,y3,j) &
      & + ARR_B(x1,x2,x3,i,y1,y2,y3,j)*ARR_A(x1,x2,x3,i,y1,y2,y3,j)
end do; end do; ...; end do

to an implementation using statement functions as abbreviation,

A(i,j) = ARR_A(x1,x2,x3,i,y1,y2,y3,j)
B(i,j) = ARR_B(x1,x2,x3,i,y1,y2,y3,j)
...
do x1=1,NGRID1; do x2=1,NGRID2; ...; do y3=1,NGRIDN
  x(i,j) = x(i,j) + A(i,j)*B(i,j) + B(i,j)*A(i,j)
end do; end do; ...; end do

In the second version it seems much easier to spot the mixed-uo order of the indices after noticing that there is an issue.

Using macros?

The only method to achieve similar semantics I found was using preprocessor macros. Typically, C-style macros are used, which here would lead to

#define A(i,j) ARR_A(x1,x2,x3,i,y1,y2,y3,j)

but this sacrifices scoping of the name. While this seems potentially useful, especially when the same abbreviation is used accross multiple functions, it seems potentially rather confusing. For a start, the compiler will give more useful error messages, when using statement functions incorrectly than it can for macros.

It also introduces the need to manually #undefine the name, if we want to reuse the same name in different contexts for different things.

Solution (1)

I had misunderstood how CONTAINS and internal functions work (see IanH's answer and its comments).

Short: Internal functions can be defined inside subroutines with CONTAINS, and work exactly like statement functions, except for a more verbose syntax.

The misunderstanding arose partly because it apparently is not allowed to define a subroutine inside a subroutine in this manner. Apparently it does work. No idea what I was doing wrong.

Answer 1

Internal procedures are pretty much a drop-in replacement for statement functions for this use case, with their usual benefits of more capability and less proneness to error, at the cost of some slight extra verbosity with their definition.

MODULE some_module     ! Just for the sake of example.
  ...
CONTAINS
  ! The host of the internal procedure.  This could 
  ! also be an external procedure or a main program.
  ! It could also be a function.
  SUBROUTINE some_subroutine
    ...
    DO i = 1, 3
      DO j = 1, 3
        DO x1 = 1, 3
          DO x2 = 1, 3
            ...
            x(i,j) = x(i,j) + A(i,j)*B(i,j) + B(i,j)*A(i,j)
            ...
          END DO
        END DO
      END DO
    END DO     
    ...      
  CONTAINS
    ! An internal procedure.  For different use cases 
    ! this could also be a subroutine.
    FUNCTION A(i,j)
      INTEGER, INTENT(IN) :: i,j
      REAL :: A
      A = ARR_A(x1,x2,x3,i,y1,y2,y3,j)
    END FUNCTION A

    FUNCTION B(i,j)
      INTEGER, INTENT(IN) :: i,j
      REAL :: B
      B = ARR_B(x1,x2,x3,i,y1,y2,y3,j)
    END FUNCTION B
    ...
 END SUBROUTINE some_subroutine
 ...
END MODULE some_module

Note that the maximum array rank under the rules of Fortran 77/90/95/2003 is seven.