What is wrong with my Fortran dynamic allocation code?

Question

This is my first question on the internet!!! Also, this is a Fortan question and I just started to learn this language. So, please forgive my ignorance. Specifically, I apologize if my example code is large. I have done my best to reduce its size to bare minimum without undermining its clarity.

Here is the problem: I am trying to create a dynamic simulation model that uses a function that will be called multiple times during a model run. The function itself is a dynamic integration system so the values of the functions must be retained at each time step of the simulation so they can be used in the next time step.

The model includes a simple stock (P) that reduces over time by a constant rate (10% per time). The function is a exponential delay. It takes 4 arguments: (input, delay time, delay order, input initial). Note that the output of the function has 3 dimensions. The first dimension represents different instances of the delay call within the program. The second dimension represents order of the delay. And the third dimension represents the simulation time.

The main challenge here was to dynamically allocate the arrays that will be created by the function. To address this issue, I have tried to follow the guidelines provided here: How to increase array size on-the-fly in Fortran? The program I have written sounds to be working fine in some cases but not always. For example, if you compile the codes as they appear below, you will get the correct result (correct behavior is that all variables should start from 100 and gradually decline). But if you change the order of delay in U from 5 to 6, then the results will be wrong (U starts from 100 but then fluctuates and sometimes go above 100).

What am I missing here? Is there anything wrong with the way I have implemented move_alloc? Or the problem is somewhere else? Could it be the compiler (gfortran) that is failing? I truly appreciate your help, in advance!

program model

  implicit none
  real, parameter :: start = 0, end = 10, dt = 0.25
  integer, parameter :: n = int((end - start)/dt)
  real, dimension(0:n+1) :: time, P, R, S, W, U
  integer :: i
  real, dimension(:,:,:), allocatable :: DN, temp

  P(0) = 100
  time(0) = start

  open (unit=10, file='out.txt', status='replace', action='write', position='rewind')
  write (10, *) 'Time, ', 'P, ', 'R, ', 'S, ', 'W, ', 'U'

  do i = 0, n
     P(i+1) = P(i) - .1*P(i)*dt
     R(i) = delay(P(i), 3.0, 4, P(0))
     S(i) = delay(P(i), 2.0, 4, P(0))
     W(i) = delay(P(i), 4.0, 5, P(0))
     U(i) = delay(P(i), 1.0, 5, P(0))
     time(i+1) = time(i) + dt
     write (10, *) time(i), ',', P(i), ',', R(i), ',', S(i), ',', W(i), ',', U(i)
  end do

  close (10)

contains
  function delay(input, dtime, order, init)
    real :: delay
    real, intent(in) :: input, dtime, init
    integer, save :: j, k
    integer, intent(in) :: order
    integer :: l
    data j /0/
    data k /0/
    if (j == i) then
       k = k + 1
       allocate(temp(1:k, 1:order, 0:n+1))
       if (allocated(DN)) temp(1:k-1, 1:order, 0:n+1) = DN
       call move_alloc(temp, DN)
    else
       k = 1
       j = i
    end if
    if (i == 0) then
       do l = 1, order
          DN(k, l, 0) = init * dtime / order
       end do
    end if
    DN(k, 1, i+1) = DN(k, 1, i) + (input - DN(k, 1, i) * order / dtime) * dt
    if (order > 1) then
       do l = 2, order
          DN(k, l, i+1) = DN(k, l, i) + (DN(k, l-1, i) - DN(k, l, i)) * order * dt / dtime
       end do
    end if
    delay = DN(k, order, i) * order / dtime
  end function delay

end program model

Answer 1

You are indexing your arrays incorrectly, and should be using compiler run-time checks to detect these errors in your program.

In the function delay you are relying on state stored in DN, and reference sections of this array using the order dummy argument. Unfortunately, order isn't tied correctly to the actual shape of DN.

Consider the lines

   allocate(temp(1:k, 1:order, 0:n+1))
   if (allocated(DN)) temp(1:k-1, 1:order, 0:n+1) = DN
   call move_alloc(temp, DN)

It is clear that after this set of lines is first executed DN is allocated and has shape [1, 4, 42]. After the second it is equally clear that its shape is [2 4 42]. Come to the third entry to the function and order is now 5 and so the shapes in the assignment do not match.

I find the code to be quite unclear so won't attempt to provide suggestions on how to correct it. In particular, the mixture of saved local state and variables from the host scope is confusing.