fortran - 从 Fortran 模块中定义的函数打印到标准输出

Question

我正在尝试学习 Fortran(不幸的是，这对我的研究小组来说是必需品)——我为自己设定的任务之一是将 Numerical Recipes 书中的一个必要函数(关联的勒让德多项式)打包到一个符合 fortran 03 的模块中。原始程序 (f77) 具有以下形式的一些错误处理:

if(m.lt.0.or.m.gt.1.or.abs(x).gt.1)pause 'bad arguments in plgndr'

自 f77 以来，Pause 似乎已被弃用，因为使用这一行会给我一个编译错误，所以我尝试了以下操作:

module sha_helper
    implicit none
    public :: plgndr, factorial!, ylm

contains
    ! numerical recipes Associated Legendre Polynomials rewritten for f03
    function plgndr(l,m,x) result(res_plgndr)
        integer, intent(in) :: l, m
        real, intent(in) :: x
        real :: res_plgndr, fact, pll, pmm, pmmp1, somx2
        integer ::  i,ll
        if (m.lt.0.or.m.gt.l.or.abs(x).gt.1) then 
            write (*, *) "bad arguments to plgndr, aborting", m, x
            res_plgndr=-10e6 !return a ridiculous value
        else
            pmm = 1.
            if (m.gt.0) then
                somx2 = sqrt((1.-x)*(1.+x))
                fact = 1.
                do i = 1, m
                    pmm = -pmm*fact*somx2
                    fact = fact+2
                end do
            end if
            if (l.eq.m) then
                res_plgndr = pmm
            else
                pmmp1 = x*(2*m+1)*pmm
                if(l.eq.m+1) then
                    res_plgndr = pmmp1
                else
                    do ll = m+2, l
                        pll = (x*(2*ll-1)*pmmp1-(ll+m-1)*pmm)/(ll-m)
                        pmm = pmmp1
                        pmmp1 = pll
                    end do
                    res_plgndr = pll
                end if
            end if
        end if
    end function plgndr

    recursive function factorial(n) result(factorial_result)
        integer, intent(in) :: n
        integer, parameter :: RegInt_K = selected_int_kind(20) !should be enough for the factorials I am using
        integer (kind = RegInt_K) :: factorial_result
        if (n <= 0) then
            factorial_result = 1
        else 
            factorial_result = n * factorial(n-1)
        end if 
    end function factorial

!     function ylm(l,m,theta,phi) result(res_ylm)
!         integer, intent(in) :: l, m
!         real, intent(in) :: theta, phi
!         real :: res_ylm, front_block
!         real, parameter :: pi = 3.1415926536
!         front_block = sqrt((2*l+1)*factorial(l-abs(m))/(4*pi*))
!     end function ylm

end module sha_helper

else 之后的主代码可以工作，但是如果我执行我的主程序并使用错误值调用该函数，则程序会在执行 print 语句之前卡住。我知道 print 语句是问题所在，因为注释掉它可以让函数正常执行，返回 -10e6 作为值。理想情况下，我希望程序在给出用户可读的错误消息后崩溃，因为给 plgndr 函数提供错误的值是程序的 fatal error 。程序 sha_lmc 正在使用函数 plgndr。目前所做的只是读取一些数组，然后打印一个 plgndr 的值进行测试(早期)。模块 sha_helper 中的函数 ylm 也没有完成，因此被注释掉了。该代码使用 gfortran sha_helper.f03 sha_lmc.f03 -o sha_lmc 编译，并且
gfortran --版本
GNU Fortran (GCC) 4.8.2

!Spherical Harmonic Bayesian Analysis testbed for Lagrangian Dynamical Monte Carlo

program sha_analysis
use sha_helper
implicit none
!Analysis Parameters
integer, parameter :: harm_order = 6
integer, parameter :: harm_array_length = (harm_order+1)**2
real, parameter :: coeff_lo = -0.1, coeff_hi = 0.1, data_err = 0.01 !for now, data_err fixed rather than heirarchical
!Monte Carlo Parameters
integer, parameter :: run = 100000, burn = 50000, thin = 100
real, parameter :: L = 1.0, e = 1.0

!Variables needed by the program
integer :: points, r, h, p, counter = 1
real, dimension(:), allocatable :: x, y, z 
real, dimension(harm_array_length) :: l_index_list, m_index_list
real, dimension(:,:), allocatable :: g_matrix

!Open the file, allocate the x,y,z arrays and read the file
open(1, file = 'Average_H_M_C_PcP_boschi_1200.xyz', status = 'old')
read(1,*) points
allocate(x(points))
allocate(y(points))
allocate(z(points))
print *, "Number of Points: ", points
readloop: do r = 1, points
    read(1,*) x(r), y(r), z(r)
end do readloop

!Set up the forwards model
allocate(g_matrix(harm_array_length,points))
!Generate the l and m values of spherical harmonics
hloop: do h = 0, harm_order
    ploop: do p = -h,h
        l_index_list(counter) = h
        m_index_list(counter) = p
        counter = counter + 1
    end do ploop
end do hloop

print *, plgndr(1,2,0.1)
!print *, ylm(1,1,0.1,0.1)


end program sha_analysis

Answer 1

您的程序执行所谓的递归 IO - 对 plgndr 的初始调用位于 IO 语句(打印语句)的输出项列表中 [将输出定向到控制台] - 在该函数内部，您还尝试执行另一个 IO 语句 [输出到控制台]。这是不允许的 - 参见 F2003 的 9.11p2 和 p3 或 F2008 的 9.12p2。

一种解决方案是将函数调用与主程序中的 io 语句分开，即

REAL :: a_temporary
...
a_temporary = plgndr(1,2,0.1)
PRINT *, a_temporary

F2008 中的其他替代方案(但不是 F2003 - 因此第一段中的 [ ] 部分)包括将函数的输出定向到不同的逻辑单元(注意 WRITE (*, ... 和 PRINT ... 引用相同的单元)。

在 F2008 中，您还可以将 WRITE 语句替换为带有消息的 STOP 语句(消息必须是常量 - 这不会让您报告有问题的值)。

无意中调用递归 IO 的可能性是某些编程风格不鼓励在函数中执行 IO 的部分原因。