How to pass many temporary arrays to a subroutine in Fortran?

Question

Suppose I'm writing a Fortran subroutine that performs one time step of a simulation, and so will be called repeatedly. To do its work, it requires many temporary arrays, which aren't meaningful or useful outside the subroutine. The user may wish to perform several simulations of different sizes during the program's execution, maybe concurrently, but the dimensions of the arrays will be the same for all the time steps of a particular simulation.

How should I handle these temporary arrays? Every way I can think of has some problem.

• If I were to use arrays local to the subroutine, they would either have to be put on the stack, which would limit how large they could be, or they would have to be allocated on the heap, which would waste time as the arrays would be allocated and deallocated every time step.

• Alternatively, the caller could allocate the arrays. However, if the caller were to pass each array as a separate argument to the subroutine, the caller has to know how many arrays there are and how big each one is. Furthermore, if the subroutine is modified in some way that changes the number or sizes of the arrays, every call to the subroutine would have to be modified.

• One technique I've seen is for the subroutine to accept a single large array and treat it as several temporary arrays concatenated together. However, this seems risky, since an arithmetic error could lead to bugs that are hard to find.

• Finally, I could put all the temporary arrays into a derived type and write a separate subroutine to allocate all the arrays and another subroutine to free them. This way, only a few subroutines have to know about the temporary arrays. However, my understanding is that, at least until recent versions of Fortran, it would only be possible for the derived type to contain pointers to the arrays. This would bring in all the disadvantages associated with using pointers (potential for bugs, aliasing, etc.).

Is there a way to do this that avoids these problems? (Or is any of these problems less of a big deal than I think?)

Answer 1

You can use the save attribute on the arrays.

subroutine step(x, t, do_allocattion)
    double precision, intent(in) :: x(:)
    double precision, intent(in) :: t
    logical, intent(in), optional :: do_allocation

    integer :: n
    double precision, allocatable, save :: work1(:), work2(:)
    n = size(x)

    if (present(do_allocation)) then
        if (do_allocation) then
            allocate(work1(n))
            allocate(work2(n))
        end if
    end if

    ! do work on x
end subroutine

EDIT: alternative solution, as I have mixed up with pointer arrays.

subroutine step(x, t)
    double precision, intent(in) :: x(:)
    double precision, intent(in) :: t

    integer :: n
    double precision, allocatable, save :: work1(:), work2(:)
    n = size(x)

    if (.not. allocated(work1)) then
        allocate(work1(n))
    end if
    if (.not. allocated(work2)) then
        allocate(work2(n))
    end if

    ! do work on x
end subroutine

EDITED: Caution: you cannot simply check for "allocatedness" as it will be undefined when the program is started and the subroutine executed for the first time. You can add another optional argument to deallocate. -> This is the case for pointer arrays.