Search code examples
fortranfortran90dynamic-allocation

Subroutine not returning correct numerical values in assumed shape array due to index renaming in the main program

The argument of my fortran 95 subroutine is an assumed shape array with intent inout:

the_subroutine(my_argument)
real, dimension(:,:), intent(inout) :: my_argument
(...)

In the main program, I have an allocatable array. I allocate it and also rename indexes. Then I call the subroutine and pass that (correctly allocated) array to the subroutine:

allocate(the_array( 5:1005 , 5:1005 ))
call the_subroutine(my_argument=the_array)

The subroutine does certain calculations and fills the array with values. In the very last line before the end of the subroutine, I check a random value:

(...)
print*, my_argument(213,126) ! I get 2.873...
end subroutine the_subroutine

Then, in the very first line after the subroutine call, I check if the value has been correctly communicated by the subroutine to the outer world, but that is not the case:

call the_subroutine(my_argument=the_array)
print*, the_array(213,126) ! I get 3.798... A completely different value.

The problem arises from having re-indexed the array in the main program as:

allocate(the_array( 5:1005 , 5:1005 ))

where max_index - min_index = 1000-1, but the subroutine "sees" the array internally as if I had declared the normal way, i.e.:

allocate(the_array( 1:1000, 1:1000))

Or simply, allocate(the_array( 1000, 1000 ))

Therefore, the element (213,126) in the internal array is in another location as in the main program array. Is there any easy way out of this?

Solution

  • Finally, I found the solution.

    First, if working in Fortran 2003 (or Fortran 95 with non-standard extensions), you may simply declare the assumed shape argument in the subroutine as ALLOCATABLE:

    subroutine the_subroutine(my_argument)
real, dimension(:,:), allocatable, intent(inout) :: my_argument

    Then the subroutine "sees" the renamed index correctly. However this is not allowed in the Fortran 95 standard.

    In Fortran 95, the most elegant way I found for this is by making use of a pointer:

    program example
implicit none
real, dimension(:,:), allocatable, target  :: the_array
real, dimension(:,:),              pointer :: the_pointer
[...]
allocate(the_array(5:1005,5:1005))
the_pointer => the_array
call the_subroutine(my_argument=the_pointer)

    And in the subroutine:

    subroutine the_subroutine(my_argument)
real, dimension(:,:), pointer :: my_argument

    Then it works perfectly. Inside the subroutine, MY_ARGUMENT is treated exactly as if it was an assumed shape array.