Ich versuche, Fortran zu lernen (das ist leider eine Notwendigkeit für meine Forschungsgruppe). Eine der Aufgaben, die ich mir gestellt habe, bestand darin, eine der erforderlichen Funktionen (Associated Legendre-Polynome) aus dem Buch Numerical Recipes in ein fortran 03-kompatibles Modul zu packen. Das ursprüngliche Programm (f77) hat eine Fehlerbehandlung in Form von Folgendem:

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

Pause scheint veraltet zu sein, da f77 bei Verwendung dieser Zeile einen Kompilierungsfehler verursacht. Daher habe ich Folgendes versucht:

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

Der Hauptcode nach "else" funktioniert, aber wenn ich mein Hauptprogramm ausführe und die Funktion mit falschen Werten aufrufe, friert das Programm ein, bevor die print-Anweisung ausgeführt wird. Ich weiß, dass die print-Anweisung das Problem ist, da durch Auskommentieren die Funktion normal ausgeführt werden kann und -10e6 als Wert zurückgegeben wird. Im Idealfall möchte ich, dass das Programm abstürzt, nachdem eine vom Benutzer lesbare Fehlermeldung ausgegeben wurde, da es ein schwerwiegender Fehler für das Programm ist, der Funktion plgndr falsche Werte zuzuweisen. Die Funktion plgndr wird vom Programm sha_lmc verwendet. Derzeit werden nur einige Arrays gelesen und anschließend ein Wert von plgndr zum Testen ausgegeben (in den frühen Tagen). Die Funktion ylm im Modul sha_helper ist ebenfalls nicht beendet, daher ist sie auskommentiert. Der Code kompiliert mit gfortran sha_helper.f03 sha_lmc.f03 -o sha_lmc und 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