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.
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