Fortran77转Fortran90中GO TO循环改写后求和异常求助

我尝试将一段Fortran 77代码转换为Fortran 90，但搞不懂原代码里的GO TO逻辑，自己改写后求和逻辑不对，无法正确累加乘积。以下是原代码和我改写的代码，求帮忙排查问题。

原Fortran 77代码 #

FUNCTION YRR(L1,L2,L,CT1,CT2,DP)

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      
      PI=dacos(-1d0)
      
      XL1=L1
      XL2=L2
      XL=L
      DEN=SQRT(4.D0*PI)*2.D0*PI
      FACT=((2.D0*XL+1.D0)/DEN)*PARITII(L1+L2)

      M=0
      
      SUM=THRJ(XL1,XL2,XL,dble(0),dble(0),dble(0))*PLM(L1,0,CT1)*PLM(L2,0,CT2)*FACT
 2000 M=M+1
      IF (M.Gt.min(l1,l2)) GO TO 3000
      SUM=SUM+2.D0*PARITii(M)*THRJ(XL1,XL2,XL,dble(M),dble(-M),dble(0)) &
              *PLM(L1,M,CT1)*PLM(L2,M,CT2)*COS(dble(M)*DP)*FACT
              
      GO TO 2000
        
 3000 YRR=SUM
      RETURN
 
      END

改写后的Fortran 90代码 #

FUNCTION YRR(L1, L2, L, CT1, CT2, DP)
  IMPLICIT NONE
  INTEGER, INTENT(IN) :: L1, L2, L
  REAL(8), INTENT(IN) :: CT1, CT2, DP
  REAL(8) :: YRR
  REAL(8), PARAMETER :: PI = ACOS(-1.0d0)
  REAL(8) :: XL1, XL2, XL, DEN, FACT
  INTEGER :: M
  REAL(8) :: SUM
  XL1 = L1
  XL2 = L2
  XL = L
  DEN = SQRT(4.0d0 * PI) * 2.0d0 * PI
  FACT = ((2.0d0 * XL + 1.0d0) / DEN) * PARITII(L1 + L2)
  M = 0
  SUM = THRJ(XL1, XL2, XL, 0.0d0, 0.0d0, 0.0d0) * PLM(L1, 0, CT1) * PLM(L2, 0, CT2) * FACT
  DO WHILE (M <= MIN(L1, L2))
    M = M + 1
    SUM = SUM + 2.0d0 * PARITii(M) * THRJ(XL1, XL2, XL, M, -M, 0.0d0) &
        * PLM(L1, M, CT1) * PLM(L2, M, CT2) * COS(M * DP) * FACT
  END DO
  YRR = SUM
  RETURN
END FUNCTION YRR

问题分析与修正 #

原代码循环逻辑拆解 #

原Fortran 77的GO TO循环实际执行流程：

1. 初始化M=0，先计算SUM的初始项（M=0的情况）
2. 进入标签2000，先将M加1
3. 判断M是否大于MIN(L1,L2)：
   • 如果是，跳转到3000结束循环
   • 如果否，执行累加操作，然后跳回2000继续循环
     核心：循环中M的取值范围是1到MIN(L1,L2)，当M等于MIN(L1,L2)时会执行累加，M加1后超过阈值才退出

改写代码的错误点 #

你用DO WHILE(M <= MIN(L1, L2))的逻辑会导致：

• 初始M=0满足条件，进入循环
• M加1变为1，执行累加
• 直到M变为MIN(L1,L2)+1时才退出循环，这意味着多执行了一次M=MIN(L1,L2)+1的循环
  此时PLM(L1, M, CT1)等函数可能因M超出合法范围返回错误值，或者累加了不该存在的项，最终导致求和结果错误。

修正后的代码 #

最简洁的方式是直接用Fortran 90的DO循环指定范围，完全匹配原代码逻辑：

FUNCTION YRR(L1, L2, L, CT1, CT2, DP)
  IMPLICIT NONE
  INTEGER, INTENT(IN) :: L1, L2, L
  REAL(8), INTENT(IN) :: CT1, CT2, DP
  REAL(8) :: YRR
  REAL(8), PARAMETER :: PI = ACOS(-1.0d0)
  REAL(8) :: XL1, XL2, XL, DEN, FACT
  INTEGER :: M
  REAL(8) :: SUM
  XL1 = L1
  XL2 = L2
  XL = L
  DEN = SQRT(4.0d0 * PI) * 2.0d0 * PI
  FACT = ((2.0d0 * XL + 1.0d0) / DEN) * PARITII(L1 + L2)
  M = 0
  SUM = THRJ(XL1, XL2, XL, 0.0d0, 0.0d0, 0.0d0) * PLM(L1, 0, CT1) * PLM(L2, 0, CT2) * FACT
  ! 直接指定M从1到MIN(L1,L2)，完全匹配原循环逻辑
  DO M = 1, MIN(L1, L2)
    SUM = SUM + 2.0d0 * PARITii(M) * THRJ(XL1, XL2, XL, DBLE(M), DBLE(-M), 0.0d0) &
        * PLM(L1, M, CT1) * PLM(L2, M, CT2) * COS(DBLE(M) * DP) * FACT
  END DO
  YRR = SUM
  RETURN
END FUNCTION YRR

如果偏好DO WHILE风格，也可以调整逻辑完全复刻原代码的跳转逻辑：

FUNCTION YRR(L1, L2, L, CT1, CT2, DP)
  IMPLICIT NONE
  INTEGER, INTENT(IN) :: L1, L2, L
  REAL(8), INTENT(IN) :: CT1, CT2, DP
  REAL(8) :: YRR
  REAL(8), PARAMETER :: PI = ACOS(-1.0d0)
  REAL(8) :: XL1, XL2, XL, DEN, FACT
  INTEGER :: M
  REAL(8) :: SUM
  XL1 = L1
  XL2 = L2
  XL = L
  DEN = SQRT(4.0d0 * PI) * 2.0d0 * PI
  FACT = ((2.0d0 * XL + 1.0d0) / DEN) * PARITII(L1 + L2)
  M = 0
  SUM = THRJ(XL1, XL2, XL, 0.0d0, 0.0d0, 0.0d0) * PLM(L1, 0, CT1) * PLM(L2, 0, CT2) * FACT
  DO
    M = M + 1
    IF (M > MIN(L1, L2)) EXIT
    SUM = SUM + 2.0d0 * PARITii(M) * THRJ(XL1, XL2, XL, DBLE(M), DBLE(-M), 0.0d0) &
        * PLM(L1, M, CT1) * PLM(L2, M, CT2) * COS(DBLE(M) * DP) * FACT
  END DO
  YRR = SUM
  RETURN
END FUNCTION YRR

内容的提问来源于stack exchange，提问作者Mekobi