type(c_ptr) Example

a.f90

program a
use types, only: dp
use compute, only: init, register_func, run, eq, destroy, get_context
use iso_c_binding, only: c_ptr, c_loc, c_f_pointer

type my_data
    ! Material coefficients:
    real(dp) :: a11, a12, a21, a22
    ! There can be a lot of variables and big arrays here, this needs
    ! to be passed around by reference.
end type

type(eq), pointer :: d
type(my_data), target :: data1, data2

data1%a11 = 0
data1%a12 = -1
data1%a21 = 1
data1%a22 = 0

data2%a11 = 0
data2%a12 = 1
data2%a21 = 1
data2%a22 = 0

call init(d)
call register_func(d, derivs, c_loc(data1))
call run(d, [0.0_dp, 1.0_dp], 0.1_dp, 10)
call print_material_parameters(d)
print *
call register_func(d, derivs, c_loc(data2))
call run(d, [0.0_dp, 1.0_dp], 0.1_dp, 10)
call print_material_parameters(d)
call destroy(d)

contains

subroutine print_material_parameters(d)
type(eq), intent(in) :: d
type(my_data), pointer :: ctx
call c_f_pointer(get_context(d), ctx)
print "('Material parameters: ', f0.6, ' ', f0.6, ' ', f0.6, ' ', f0.6)", &
    ctx%a11, ctx%a12, ctx%a21, ctx%a22
end subroutine

function derivs(x, data) result(y)
use types, only: dp
real(dp), intent(in) :: x(2)
type(c_ptr), intent(in) :: data
real(dp) :: y(2)
type(my_data), pointer :: d
call c_f_pointer(data, d)
y(1) = d%a11 * x(1) + d%a12 * x(2)
y(2) = d%a21 * x(1) + d%a22 * x(2)
end function

end program

compile

#! /bin/bash

set -e
FFLAGS="-Wall -Wextra -Wimplicit-interface -fPIC -Werror -fmax-errors=1 -g -fbounds-check -fcheck-array-temporaries -fbacktrace"
gfortran $FFLAGS -c types.f90 -o types.o
gfortran $FFLAGS -c compute.f90 -o compute.o
gfortran $FFLAGS -c a.f90 -o a.o
gfortran $FFLAGS -o a a.o compute.o types.o

compute.f90

module compute
use types, only: dp
use iso_c_binding, only: c_ptr
implicit none
private
public init, destroy, register_func, run, eq, get_context

abstract interface
    function derivs(x, data) result(y)
    use types, only: dp
    use iso_c_binding, only: c_ptr
    implicit none
    real(dp), intent(in) :: x(2)
    type(c_ptr), intent(in) :: data
    real(dp) :: y(2)
    end function
end interface

type eq
    type(c_ptr) :: data
    procedure(derivs), nopass, pointer :: func
end type

contains

subroutine init(d)
type(eq), pointer, intent(inout) :: d
allocate(d)
d%func => NULL()
end subroutine

subroutine destroy(d)
type(eq), pointer, intent(inout) :: d
deallocate(d)
end subroutine

subroutine register_func(d, func, data)
type(eq), intent(inout) :: d
procedure(derivs) :: func
type(c_ptr), intent(in) :: data
d%func => func
d%data = data
end subroutine

function get_context(d) result(data)
type(eq), intent(in) :: d
type(c_ptr) :: data
data = d%data
end function

subroutine run(d, x0, dt, n_steps)
type(eq), intent(in) :: d
real(dp), intent(in) :: x0(2), dt
integer, intent(in) :: n_steps

real(dp) :: x(2), dx(2), t
integer :: i
if (.not. associated(d%func)) then
    print *, "d%func is not associated"
end if
x = x0
t = 0
do i = 1, n_steps
    dx = d%func(x, d%data)
    print "(f10.6, f10.6)", x
    x = x + dx * dt
    t = t + dt
end do
end subroutine

end module

types.f90

module types

implicit none
private
public dp

integer, parameter :: dp=kind(0.d0)

end module