Long time ago I designed a generic ODE solver for the problem:
dy/dt = f(y,t) y(t0)=y0 t=[t0,tmax]
Until now my design (see below) has worked OK, but now I have
encountered a funny problem. I want to include the NETLIB package
'rksuite_90' into my ode_handler. This package has a slightly
different argument list for the function f. ( f(t,y) instead of f(y,t). )
Is it possible to include 'rksuite_90' into my module procedure
'ode_integrate' without modifying the argument lists for f in my program
or the NETLIB routines? (Modification of argument lists is impossible
if I had to include two different commercial libraries.)
thanks for any help,
Jens Helmers
A simplified version of my "generic" integrator is:
module ode_handler
....
contains
....
subroutine ode_integrate( method, t0, tmax, y, f )
character(len=*), intent(in) :: method
real, intent(in) :: t0, tmax
real, dimension(:), intent(inout) :: y ! y0 on input, y(tmax) on exit
interface
function f(y,t)
real, dimension(:), intent(in) :: y
real, intent(in) :: t
real, dimension(size(y)) :: f
end function fun
end interface
select case(method)
case('euler_explicit')
call euler_explicit(f,t0,tmax,y)
case('euler_implicit')
call euler_implicit(f,t0,tmax,y)
case('rksuite_90')
call rk_suite_90(f,t0,tmax,y) ! Interface error for f
call rk_suite_90(f_alias,t0,tmax,y) ! Not legal in f95,f2k
case('...')
...
end select
...
contains
function f_alias(t,y)
real, dimension(:), intent(in) :: y
real, intent(in) :: t
real, dimension(size(y)) :: f_alias
f_alias = f(y,t)
end function f_alias
end subroutine ode_integrate
end module ode_handler
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|