Why a function returns an array is much slower than a subroutine returns an array? (real MWE included)

[shahmoradi]

That is likely because you have not yet hit the stack limit. GFortran has -fmax-stack-var-size=n to request allocations on the heap for arrays of size n bytes or larger.