Fortran动态分配代码有什么问题？

Question

这是我在互联网上的第一个问题！！！另外，这是一个Fortan问题，我才刚开始学习这种语言。所以，请原谅我的无知。具体来说，如果示例代码过大，我深表歉意。我已尽力将其尺寸减小到最小，而不会破坏其清晰度。

这里是问题：我正在尝试创建一个动态仿真模型，该模型使用在模型运行期间会多次调用的函数。该函数本身是一个动态集成系统，因此必须在模拟的每个时间步都保留这些函数的值，以便可以在下一时间步中使用它们。

该模型包括一个简单的库存（P），该库存会随着时间的推移以恒定的速率减少（每次10％）。该函数是指数延迟。它包含4个参数：（输入，延迟时间，延迟顺序，输入初始值）。请注意，函数的输出具有3个维度。第一维表示程序内延迟调用的不同实例。第二维表示延迟的顺序。第三维代表模拟时间。

这里的主要挑战是动态分配将由该函数创建的数组。为了解决此问题，我尝试遵循此处提供的指导：How to increase array size on-the-fly in Fortran?我编写的程序在某些情况下听起来不错，但并非总是如此。例如，如果按如下所示编译代码，则将获得正确的结果（正确的行为是所有变量应从100开始并逐渐减小）。但是，如果将U的延迟顺序从5更改为6，则结果将是错误的（U从100开始，然后波动，有时会超过100）。

我在这里想念什么？我实现move_alloc的方式有什么问题吗？还是问题出在别的地方？可能是编译器（gfortran）失败了吗？提前感谢您的帮助！

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

在函数delay中，您依赖于存储在DN中的状态，并使用order虚拟参数引用此数组的部分。不幸的是，order没有正确绑定到DN的实际形状。

考虑线

   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)

很明显，在第一次执行这组线之后，DN被分配并具有形状[1, 4, 42]。在第二秒之后，同样清楚的是其形状为[2 4 42]。转到函数的第三个条目，order现在为5，因此赋值中的形状不匹配。

我发现代码尚不清楚，因此不会尝试提供有关纠正方法的建议。特别是，主机作用域中保存的局部状态和变量的混合令人困惑。