Passing an internal procedure as argument

Question

I want to solve a differential equation lots of times for different parameters. It is more complicated than this, but for the sake of clarity let's say the ODE is y'(x) = (y+a)*x with y(0) = 0 and I want y(1). I picked the dverk algorithm from netlib for solving the ODE and it expects the function on the right hand side to be of a certain form. Now what I did with the Intel Fortran compiler is the following (simplified):

subroutine f(x,a,ans)
implicite none
double precision f,a,ans,y,tol,c(24),w(9)
...
call dverk(1,faux,x,y,1.d0,tol,ind,c,1,w)
...
contains
    subroutine faux(n,xx,yy,yprime)
           implicite none
           integer n
           double precision xx,yy(n),yprime(n)
           yprime(1) = (yy(1)+a)*xx
    end subroutine faux
end subroutine f

This works just fine with ifort, the sub-subroutine faux sees the parameter a and everything works as expected. But I'd like the code to be compatible with gfortran, and with this compiler I get the following error message:

Error: Internal procedure 'faux' is not allowed as an actual argument at (1)

I need to have the faux routine inside f, or else I don't know how to tell it the value of a, because I can't change the list of parameters, since this is what the dverk routine expects.

I would like to keep the dverk routine and understand how to solve this specific problem without a workaround, since I feel it will become important again when I need to integrate a parameterized function with different integrators.

Answer 1

You could put this all in a module, and make a a global module variable. Make faux a module procedure. That way, it has access to a.

module ode_module
    double precision::a

    contains

    subroutine f(x,a,ans)
        implicit none
        double precision f,ans,y,tol,c(24),w(9)

        call dverk(1,faux,x,y,1.d0,tol,ind,c,1,w)

    end subroutine 

    subroutine faux(n,xx,yy,yprime)
       implicite none
       integer n
       double precision xx,yy(n),yprime(n)
       yprime(1) = (yy(1)+a)*xx
    end subroutine faux

end module