title: Using an external function in a module
body: <p>I am writing a n-dimensional numerical solver in Fortran. I created a module for the same which is called from the main program. I used <code>external</code> for calling the unknown function when I wrote for first order ode. However on replicating the result for more than one dimension using <code>external</code> I get the following error </p>
  
  <pre><code>Error: EXTERNAL attribute conflicts with DIMENSION attribute in 'f_x' 
  </code></pre>
  
  <p>This means my compiler cannot process an external function with array output. </p>
  
  <p>I tried using an interface block to predefine the function parameters with fixed dimensions and I end up with this error </p>
  
  <pre><code>Error: PROCEDURE attribute conflicts with INTENT attribute in 'f_x' 
  </code></pre>
  
  <p>Is there any way in Fortran to have an external function returning an array being initialized using a dummy function in a subroutine.</p>
  
  <p>Here is the code that I used.</p>
  
  <pre><code>module int_adaptive
  implicit none
  contains
  
  subroutine RK4nd(f_x, norder, nvar, x_0, t_0, t_f, x_out, t_out)
      INTERFACE
          FUNCTION f_x (x_0, t_0)
              integer, parameter  :: dp=kind(0.d0)
              REAL(dp), INTENT(IN)    :: x_0(2, 3), t_0
              REAL, intent(out)   :: f_x(2,3)
          END FUNCTION f_x
      END INTERFACE
  
  
      integer, parameter  :: dp=kind(0.d0)
      integer             :: norder, i, nvar
      external            :: f_x
      real(dp), intent(in):: x_0(norder, nvar), t_0, t_f
      real(dp)            :: x_out(norder, nvar), t_out, k1(norder, nvar), y1(norder, nvar), y(norder, nvar)
      real(dp)            :: k2(norder, nvar), k3(norder, nvar), k4(norder, nvar),&amp;
                             y2(norder, nvar), y3(norder, nvar), y4(norder, nvar)!, f_x(norder, nvar)
      real(dp)            :: h_min, h_0, h_max, tol, err_est, h, t
  
      if (h_0&lt;h) then
          h = h_0
      end if
      if ((t_f - t_0) &lt; 0.0) then
          h = -h
      end if
      t = t_0
      y = x_0
  
      do while (t .NE. t_f)
          k1 = f_x(y, t)
          y1 = y + k1*h/2.0
  
          k2 = f_x(y1, t+h/2.0)
          y2 = y + k2*h/2.0
  
          k3 = f_x(y2, t+h/2.0)
          y3 = y + k3*h
  
          k4 = f_x(y3, t+h)
          y4 = y + (k1+ 2*k2 + 2*k3 +k4)*h/6.0
  
          t = t + h
          y = y4
  
      end do
  
      x_out = y
      t_out = t
  
  end subroutine
  end module
  </code></pre>
  
  <p>I could define a standard integrator with fixed function input inside a module and call it by name directly, but since I deal with ODE's of various orders it will be complicated each time I am changing the order I'll need to update this module as well.</p>