Print to standard output from a function defined in an Fortran module

Question

I'm trying to learn Fortran (unfortunately a necessity for my research group) - one of the tasks I set myself was to package one of the necessary functions (Associated Legendre polynomials) from the Numerical Recipes book into a fortran 03 compliant module. The original program (f77) has some error handling in the form of the following:

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

Pause seems to have been deprecated since f77 as using this line gives me a compiling error, so I tried the following:

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

The main code after the else works, but if I execute my main program and call the function with bad values, the program freezes before executing the print statement. I know that the print statement is the problem, as commenting it out allows the function to execute normally, returning -10e6 as the value. Ideally, I would like the program to crash after giving a user readable error message, as giving bad values to the plgndr function is a fatal error for the program. The function plgndr is being used by the program sha_lmc. Currently all this does is read some arrays and then print a value of plgndr for testing (early days). The function ylm in the module sha_helper is also not finished, hence it is commented out. The code compiles using gfortran sha_helper.f03 sha_lmc.f03 -o sha_lmc, and gfortran --version 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

Your program does what is known as recursive IO - the initial call to plgndr is in the output item list of an IO statement (a print statement) [directing output to the console] - inside that function you then also attempt to execute another IO statement [that outputs to the console]. This is not permitted - see 9.11p2 and p3 of F2003 or 9.12p2 of F2008.

A solution is to separate the function invocation from the io statement in the main program, i.e.

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

Other alternatives in F2008 (but not F2003 - hence the [ ] parts in the first paragraph) include directing the output from the function to a different logical unit (note that WRITE (*, ... and PRINT ... reference the same unit).

In F2008 you could also replace the WRITE statement with a STOP statement with a message (the message must be a constant - which wouldn't let you report the problematic values).

The potential for inadvertently invoking recursive IO is part of the reason that some programming styles discourage conducting IO in functions.

Answer 2

Try:

if (m.lt.0.or.m.gt.l.or.abs(x).gt.1) then
   write (*, *) "bad arguments to plgndr, aborting", m, x
   stop
else
   ...
end if