Fortran 数组自动更改值

Question

我正在做关于在 ANSYS CFX 上进行风力涡轮机模拟的硕士论文。要计算参数，我必须使用 FORTRAN 编写的代码。我遇到了一个非常奇怪的问题，需要你的帮助。完整的 fortran 代码很长，但我只会在导致问题的代码中发布子例程。

我正在使用子程序将叶片划分为径向元素。使用两种类型的划分，然后将元素的径向位置存储在两个一维数组 RI 和 RJ 中。

稍后在代码中，我遇到了一些错误，我将其追溯到这个子程序，它没有给我正确的 RJ 值。子程序给出了正确的 RI 值。然后我使用计算correct values of RJ 的相同循环显示了 RJ 的结果。幸运的是，RJ 的值是正确的。然后在同一个子程序中，在那个循环之后，我立即开始另一个循环以再次显示 RJ 的值，不幸的是这次它们不正确incorrect values of RJ，即使两个循环之间没有任何东西可以改变 RJ 的值。似乎 RJ 的每个元素都被 RI 的下一个元素替换。因此，我开始知道这是错误的根源。我重新检查了程序，似乎没有错误，错误源于子程序内部的某个地方。我在没有 ANSYS CFX 的 fortran 上使用另一个程序使用相同的子程序来计算风力涡轮机的性能，尽管这两个程序中的两个子程序完全相同，但这个问题并没有出现在该程序中。我在这件事上需要帮助，因为过去几天我都陷入了困境。我正在使用 INTEL COMPOSER 附带的 Intel fortran 编译器。

子程序如下

CALL INIT(M,PI,PREC,R,HUBRAD,RI,RJ)

SUBROUTINE INIT(M,PI,PREC,R,HUBRAD,RI,RJ)


  INTEGER:: M
  REAL:: PI,PREC,R,HUBRAD
  REAL:: RI(41),RJ(41)
  !,RRI,RRJ)

    CALL LLPOINTS (PI,PREC,R,HUBRAD,M,RI,RJ)


    return
    End Subroutine

    SUBROUTINE LLPOINTS (PI,PREC,R,HUBRAD,M,RI,RJ)

  !Input arguments: M,Pi, PREC, R, HUBRAD 
  !Output arguments : RI, RJ

  INTEGER:: M
  REAL::PREC,R,HUBRAD
  REAL:: RI(41),RJ(41)
  INTEGER :: J
  character*100 ::string1, string2, string3

  CALL MESAGE( 'WRITE', 'subroutine INIT START' )  

  CALL MESAGE( 'WRITE', "RI   RJ    DJ")

  DO J=1,M+1   


        IF (J.LT.M+1) THEN   
        RI(J)=0.5E+00*(1.E+00+HUBRAD/R)-0.5E+00*(1.E+00-HUBRAD/R)
       &*COS((J-0.5E+00)*PI/M)

        RJ(J)=0.5E+00*(1.E+00+HUBRAD/R)-0.5E+00*(1.E+00-HUBRAD/R)
       &*COS((J-1.E+00)*PI/M)


        IF (ABS(RI(J)).LT.(1/PREC)) THEN

              STOP 
        ENDIF

        ELSE
           RJ(J)=1.E+00

          END IF


          write (string1,*) RI(J)
          write (string2,*) RJ(J)


          CALL MESAGE( 'WRITE', string1//' '//string2) 

      END DO

      CALL MESAGE( 'WRITE','ri     rj')
      do j =1,m

          write (string1,*) RI(J)
          write (string2,*) RJ(J)
          CALL MESAGE( 'WRITE', string1//' '//string2 ) 

      end do
      CALL MESAGE( 'WRITE', 'subroutine LLINE OK' )  

  END SUBROUTINE

您会想知道为什么 LLPOINTS 无缘无故地存在于 INIT 子例程中。实际上在原程序中子程序INIT还有很多其他的子程序。我只使用了 LLPOINTS 并丢弃了其余的因为我不需要它们

我已经检查过我的程序在任何子程序中都没有隐含的 none，根据我的有限知识，这些子程序应该存在。当我在所有子例程中添加隐式 none 时，出现了一些编译错误，并且似乎在这些错误中，我发现存在与许多变量的声明有关的问题。我已经解决了这些问题，现在当所有问题都得到解决时，我在编译中遇到了以下错误，这在implicit none之前没有出现。

错误：函数 ACD_Dp 中引用的外部符号 LLINE 无法解析。

我不知道如何处理这个错误

Answer 1

奇怪的行为可能源于一些（不正确的）实际参数，而不是子程序本身。为了解释这一点，我们首先考虑INIT()和LLPOINTS()的简化版本如下：

subroutine INIT ( a, b )
    implicit none
    real :: a( 5 ), b( 5 )

    call LLPOINTS ( a, b )
end subroutine

subroutine LLPOINTS ( a, b )
    implicit none
    real :: a( 5 ), b( 5 )
    integer :: i

    print *, "output (1):"
    do i = 1, 5
        a( i ) = i                !! set some values to a(:) and b(:)
        b( i ) = i * 100
        print *, a( i ), b( i )   !! check the values
    enddo

    print *, "output (2):"
    do i = 1, 4
        print *, a( i ), b( i )   !! check the values again
    enddo
end subroutine

该程序将一些值设置为a(:) 和b(:) 并打印两次它们的值以进行双重检查（如在 OP 的程序中）。现在我们考虑主程序：

program main
    real :: a( 5 ), b( 5 )
    call INIT ( a, b )
end

给出了预期的结果（ifort test.f90 和 v14.0）：

 output (1):
   1.000000       100.0000
   2.000000       200.0000
   3.000000       300.0000
   4.000000       400.0000
   5.000000       500.0000
 output (2):
   1.000000       100.0000
   2.000000       200.0000
   3.000000       300.0000
   4.000000       400.0000

接下来，让我们假设 a 和 b 被错误地声明为标量变量

program main
    real :: a, b
    call INIT ( a, b )
end

甚至没有声明（即使用默认的implicit real(a-h,o-z) 规则）

program main
    call INIT ( a, b )
end

然后我们得到

 output (1):
   1.000000       100.0000
   2.000000       200.0000
   3.000000       300.0000
   4.000000       400.0000
   5.000000       500.0000
 output (2):
   1.000000       2.000000
   2.000000       3.000000
   3.000000       4.000000
   4.000000       5.000000

其模式似乎与 OP 的输出非常相似（即，输出 (2) 中的所有元素都移动了 1）。也就是说，这种奇怪行为的原因可能是我们传递了标量并访问了无效的内存区域（假设a 和b 在内存中连续对齐，并带有一些尾随内存区域）。如果是这样，main() 和 LLPOINTS() 之间的内存映射可能如下所示：

a       b      NG     NG     NG     NG
---------------------------------------
a(1)    a(2)   a(3)   a(4)   a(5)
        b(1)   b(2)   b(3)   b(4)   b(5)

如果是这种情况，由于a(i+1) = b(i)，我们在上面的输出（2）中得到了奇怪的结果。我们可以通过插入类似的行来确认这一点

if ( loc( a(2) ) == loc( b(1) ) ) stop "trapped (ifort)"

在LLPOINTS()。更重要的是，如果我们附上

ifort -check test.f90

选项，我们可以自动检测到这一点（有分段错误）。所以你能试试这个选项，看看是否是这种情况......？

Answer 2

是的，你是对的。我已经放了以下几行

If ( loc( RJ(1) ) == loc( RI(2) ) ) then
call mesage ('write','trapped (ifort)')

stop

End If

我得到了错误，打印被困 (ifort)。 但问题是，我没有在任何地方声明标量参数 RI、RJ。我一直将它们声明为向量 RI(41) 和 RJ(41)