title: Print to standard output from a function defined in an Fortran module
body: <p>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:</p>
  
  <pre><code>if(m.lt.0.or.m.gt.1.or.abs(x).gt.1)pause 'bad arguments in plgndr'
  </code></pre>
  
  <p>Pause seems to have been deprecated since f77 as using this line gives me a compiling error, so I tried the following: </p>
  
  <pre><code>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 &lt;= 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
  </code></pre>
  
  <p>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  </p>
  
  <pre><code>!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
  </code></pre>