MUSCLE and Fortran

This section shows how to use MUSCLE3 from Fortran, based on the same reaction-diffusion model given in the Python tutorial. The source code for the examples in this section is here. You can also go to docs/source/examples/fortran in the source directory (that’s easiest, and has some handy scripts and a Makefile), or copy-paste the code from here. If you’re cloning from GitHub, be sure to use the master branch, as develop may have changes incompatible with the current release.

Building and running the examples

If you’ve just built and installed the C++ version of libmuscle (which includes the Fortran bindings), then you’re all set to build the examples. To do that, go into the docs/source/examples subdirectory, activate the installation, and run Make:

~/muscle3_source$ cd docs/source/examples
~/muscle3_source$ . /path/to/muscle3/bin/muscle3.env
~/muscle3_source/docs/source/examples$ make fortran

We also need the Python version of MUSCLE3 installed, because the MUSCLE Manager comes with that, and we need it to run the simulation. The command above will also set up a Python virtual env with everything needed.

You can then run the examples using the provided scripts in docs/source/examples:

. python/build/venv/bin/activate
muscle_manager --start-all rd_implementations.ymmsl rd_fortran.ymmsl rd_settings.ymmsl
muscle_manager --start-all rd_implementations.ymmsl rd_python_fortran.ymmsl rd_settings.ymmsl

Log output

As for the Python model, the manager will make a run directory into which log files and output are written. You can set muscle_remote_log_level and muscle_local_log_level for Fortran models as you can for Python, but in Fortran this will only affect log messages generated by libmuscle, not anything written by your model. Of course, you can read muscle_local_log_level or a custom setting yourself in your F_INIT, and set your log library’s level accordingly.

If the model code logs to standard output or standard error, then the messages will end up in <rundir>/instances/<instance-id>/stdout.txt or stderr.txt. If the model logs to a file in the working directory, then you’ll find it in <rundir>/instances/<instance-id>/workdir/.

Describing Fortran implementations

The only difference between this and the Python-only example is that the Fortran implementations are used. These differ a bit from the Python versions, but they are the same as the C++ ones:

reaction_cpp:
  env:
    +LD_LIBRARY_PATH: :<muscle3_prefix>/lib
  executable: <muscle3_src>/docs/source/examples/fortran/build/reaction

diffusion_cpp:
  env:
    +LD_LIBRARY_PATH: :<muscle3_prefix>/lib
  executable: <muscle3_src>/docs/source/examples/fortran/build/diffusion

Note that <muscle3_prefix> and <muscle3_src> are not in there literally, they are just a placeholder here. When you ran Make, it created this file and substituted in your local paths. For a real model, you should put in the absolute path to the installation yourself.

Since we have compiled binaries here which can be executed directly, we don’t need to run Python, we can just specify their location directly. We do need to make sure that the MUSCLE3 libraries, which the executables are linked against, can be found. We do that by setting the LD_LIBRARY_PATH to point to their location. Since LD_LIBRARY_PATH is a list, and overwriting it may break something, we append to it by adding a + in front of the variable name. Note the initial colon in the value: this is the separator used in LD_LIBRARY_PATH, and MUSCLE3 will not add it automatically, so it needs to be there.

Reaction-diffusion in Fortran

This page assumes that you’ve at least read the Python tutorial; we’ll do the same thing again but now in Fortran.

The Fortran version of libmuscle doesn’t support running multiple components in a single program, so we’re going to do the distributed version of the reaction-diffusion example in Fortran. Here’s the basic reaction model, without MUSCLE support, in Fortran:

Reaction model

A simple reaction model on a 1D grid, in Fortran.

program reaction
    implicit none

    real (selected_real_kind(15)) :: t_cur, t_max, dt, x_max, dx, k
    integer :: i, U_size
    real (selected_real_kind(15)), dimension(:), allocatable :: U


    ! F_INIT
    t_max = 2.469136e-07
    dt = 2.469136e-08
    x_max = 1.0
    dx = 0.01
    k = -40500.0

    U_size = int(x_max / dx)
    allocate (U(U_size))
    U = 1e-20
    U(26) = 2.0
    U(51) = 2.0
    U(76) = 2.0

    t_cur = 0.0
    do while (t_cur + dt < t_max)
        ! O_I

        ! S
        U = k * U * dt
        t_cur = t_cur + dt
    end do

    ! O_F

    deallocate (U)

end program reaction

This is one of the simplest possible computational models: we set some parameters, create a state variable U, initialise the state, and then update the state in a loop until we reach the final simulation time.

The MUSCLE version in Fortran looks quite similar to the Python version:

docs/source/examples/fortran/reaction.f90

program reaction
    use ymmsl
    use libmuscle
    implicit none

    type(LIBMUSCLE_PortsDescription) :: ports
    type(LIBMUSCLE_Instance) :: instance

    type(LIBMUSCLE_Message) :: rmsg
    type(LIBMUSCLE_DataConstRef) :: rdata, item


    type(LIBMUSCLE_Message) :: smsg
    type(LIBMUSCLE_Data) :: sdata

    real (selected_real_kind(15)) :: t_cur, t_max, dt, k
    real (selected_real_kind(15)), dimension(:), allocatable :: U


    ports = LIBMUSCLE_PortsDescription()
    call ports%add(YMMSL_Operator_F_INIT, 'initial_state')
    call ports%add(YMMSL_Operator_O_F, 'final_state')
    instance = LIBMUSCLE_Instance(ports)
    call LIBMUSCLE_PortsDescription_free(ports)

    do while (instance%reuse_instance())
        ! F_INIT
        t_max = instance%get_setting_as_real8('t_max')
        dt = instance%get_setting_as_real8('dt')
        k = instance%get_setting_as_real8('k')

        rmsg = instance%receive('initial_state')
        rdata = rmsg%get_data()
        allocate (U(rdata%size()))
        call rdata%elements(U)
        call LIBMUSCLE_DataConstRef_free(rdata)

        t_cur = rmsg%timestamp()
        t_max = rmsg%timestamp() + t_max
        call LIBMUSCLE_Message_free(rmsg)

        do while (t_cur + dt < t_max)
            ! O_I

            ! S
            U = k * U * dt
            t_cur = t_cur + dt
        end do

        ! O_F
        sdata = LIBMUSCLE_Data_create_grid(U, 'x')
        smsg = LIBMUSCLE_Message(t_cur, sdata)
        call instance%send('final_state', smsg)
        call LIBMUSCLE_Message_free(smsg)
        call LIBMUSCLE_Data_free(sdata)
        deallocate (U)
    end do

    call LIBMUSCLE_Instance_free(instance)

end program reaction

We’ll go through it top to bottom, one piece at a time.

Modules

use ymmsl
use libmuscle
implicit none

Here we tell Fortran that we’ll be using the ymmsl and libmuscle modules. These mirror the corresponding Python packages. Like in C++, yMMSL support is limited to what is needed to implement components. Loading, manipulating and saving yMMSL documents is better done in Python.

Variables

type(LIBMUSCLE_PortsDescription) :: ports
type(LIBMUSCLE_Instance) :: instance

type(LIBMUSCLE_Message) :: rmsg
type(LIBMUSCLE_DataConstRef) :: rdata, item


type(LIBMUSCLE_Message) :: smsg
type(LIBMUSCLE_Data) :: sdata

real (selected_real_kind(15)) :: t_cur, t_max, dt, k
real (selected_real_kind(15)), dimension(:), allocatable :: U

Next, it’s time to declare our variables. We have a few extra variables here compared to the non-MUSCLE version. The LIBMUSCLE_PortsDescription is used to describe the ports we’ll use to send and receive to MUSCLE3. We need a LIBMUSCLE_Instance as the main interface to MUSCLE. We’ll receive messages into rmsg, a LIBMUSCLE_Message variable, and extract the data from them into rdata and item. Since received data cannot be modified, this variable is of type LIBMUSCLE_DataConstRef. For sending messages we have smsg, and we’ll create a LIBMUSCLE_Data object to put our data into before sending. Note that the names are all prefixed with LIBMUSCLE, to make sure that they don’t collide with any other library you may want to use.

Eagle-eyed readers will have noticed that dx and x_max are missing. That is because we’ll derive the size of the state vector of the model (U) from the state we receive, rather than from the configuration.

Creating an Instance

ports = LIBMUSCLE_PortsDescription()
call ports%add(YMMSL_Operator_F_INIT, 'initial_state')
call ports%add(YMMSL_Operator_O_F, 'final_state')
instance = LIBMUSCLE_Instance(ports)
call LIBMUSCLE_PortsDescription_free(ports)

In order to talk to the rest of the MUSCLE simulation, we need a LIBMUSCLE_Instance object, and we need to pass a description of the ports we’ll use when we make that. So here, we first create a LIBMUSCLE_PortsDescription, and add some ports to the F_INIT and O_F operators. We can then create the Instance object, passing the ports. Finally, we need to free the PortsDescription object. You will have to be careful to always free any object you create, or that is returned by a MUSCLE function, after you’re done using it. So we free the PortsDescription here, but we don’t free the Instance, since we need it in the following part of the code.

Note

Although the Fortran 2003 standard defines finalizers, these are not fully implemented by all compilers (notably gfortran has an incomplete implementation, see this bug tracker). Call the free() functions to ensure your program does not leak memory.

Reuse loop and settings

do while (instance%reuse_instance())
    ! F_INIT
    t_max = instance%get_setting_as_real8('t_max')
    dt = instance%get_setting_as_real8('dt')
    k = instance%get_setting_as_real8('k')

As in Python, we have a reuse loop, and except for the syntax differences between the languages, the code is exactly the same. Getting settings works as shown; with Fortran being a statically typed language, we need to specify the type we’re expecting by calling a different function. If the actual type is different, an error will be printed and the program will be halted. It is possible to detect and handle errors instead, see below for that.

Getting settings is done via the instance%get_setting_as_<type>() functions. Supported types are character, logical, int8 (a 64-bit integer), real8 (a 64-bit double precision real number), real8array (a 1D array of double precision reals) and real8array2 (a 2D array of double precision reals). See the API documentation for details.

Receiving messages and DataConstRef

Instead of initialising our state from a constant, we’re going to receive a message on our F_INIT port with the initial state:

rmsg = instance%receive('initial_state')
rdata = rmsg%get_data()
allocate (U(rdata%size()))
call rdata%elements(U)
call LIBMUSCLE_DataConstRef_free(rdata)

Calling the LIBMUSCLE_Instance_receive() function with an instance and a port name yields an object of type LIBMUSCLE_Message containing the received message. As always when using MUSCLE from Fortran, we have to remember to free the returned Message object when we are done with it. That’s not yet the case though, because we still need to get the received data out. In Python, Message is a very simple class with several public members. In Fortran, objects are always accessed and manipulated through LIBMUSCLE functions, in this case LIBMUSCLE_Message_get_data(), which we use to get a DataConstRef object. Again, since MUSCLE gave us an object, we have to remember to free it when we’re done.

First though, we’ll use the data to initialise our state. We are going to assume that we’ll receive a 1D grid of numbers, just like in the equivalent examples in the other supported languages. LIBMUSCLE_DataConstRef_size() will return the size of the array that we received (or print an error and stop if we received something else that doesn’t have a size, see below for handling errors differently). We use that to allocate the U array to the same size as the received state, and then we copy the elements of the array into U using LIBMUSCLE_DataConstRef_elements(). We can then free the rdata object, since we don’t need it any longer. We’ll hold on to the rmsg a bit longer, since we’ll need it later.

If all this freeing objects is starting to sound tedious, that’s because it is, and it’s why more modern languages like Python and C++ do this for you. Unfortunately, for Fortran, it has to be done manually.

Note that indices for the received array start at 1, as usual in Fortran. MUSCLE 3 follows the language in which you’re using it and automatically translates, so if this grid was sent from Python or C++, then received item 1 corresponds to sent item 0.

If you have a DataConstRef object, then you can check which kind of value it contains, e.g. using LIBMUSCLE_DataConstRef_is_a_grid_of_real8(). Here, we don’t bother with a check. Instead, we blindly assume that we’ve been sent a 1D grid of doubles. If that’s not the case, an error will be printed and our program will halt. That’s okay, because it means that there’s something wrong somewhere that we need to fix. MUSCLE is designed to let you get away with being a bit sloppy as long as things actually go right, but it will check for problems and let you know if something goes wrong. If you want to make a submodel or component that can handle different kinds of messages, then these inspection functions will help you do so however.

t_cur = rmsg%timestamp()
t_max = rmsg%timestamp() + t_max
call LIBMUSCLE_Message_free(rmsg)

do while (t_cur + dt < t_max)
    ! O_I

    ! S
    U = k * U * dt
    t_cur = t_cur + dt
end do

The main loop of the model is almost unchanged, we just get the starting time from the received message’s timestamp so as to synchronise with the overall simulation time correctly. The stopping time is adjusted accordingly as well. Note that t_cur and t_max are numbers, not objects, so they do not have to be freed. rmsg however does, and since we have all the information we need from the received message, we free it here.

Sending messages and Data

    ! O_F
    sdata = LIBMUSCLE_Data_create_grid(U, 'x')
    smsg = LIBMUSCLE_Message(t_cur, sdata)
    call instance%send('final_state', smsg)
    call LIBMUSCLE_Message_free(smsg)
    call LIBMUSCLE_Data_free(sdata)
    deallocate (U)
end do

call LIBMUSCLE_Instance_free(instance)

Having computed our final state, we will send it to the outside world on the final_state port. In this case, we need to send a vector of doubles, which we first need to wrap up into a Data object. A Data object works just like a DataConstRef, except that it isn’t constant, and can thus be modified. We do this by creating a Data object containing a grid value, with the array U and the index name x.

Having put our data into the Data object is a grid, we can then put the grid into a new LIBMUSCLE_Message object, and call the LIBMUSCLE_Instance_send() function to send it. Finally, we free the Message and Data objects that we created, and deallocate the state as in the original non-MUSCLE version.

That concludes the reuse loop. When we’re done running the model, the reuse loop will finish, and we can free our Instance object before we quit.

LIBMUSCLE_Data is quite versatile, and makes it easier to send data of various types between submodels. Here are some other examples of creating LIBMUSCLE_Data objects containing different kinds of data (note that freeing objects is omitted here for brevity, of course you have to do that in your model!):

type(LIBMUSCLE_Data) :: d1, d2, d3, d4, d5, d6, d7, d8, d9, d10, d11, d12
integer (LIBMUSCLE_int8), dimension(2, 3) :: ar
character(len=1), dimension(1024) :: bytes

! create a Data containing a nil value (None in Python)
d1 = LIBMUSCLE_Data()

! character strings, logicals
d2 = LIBMUSCLE_Data('String data')
d3 = LIBMUSCLE_Data(.true.)

! various kinds of numbers
d4 = LIBMUSCLE_Data(0)                   ! integer
d5 = LIBMUSCLE_Data(100_LIBMUSCLE_int8)  ! 64-bit integer
d6 = LIBMUSCLE_Data(3.141592_LIBMUSCLE_real4)    ! 32-bit real
d7 = LIBMUSCLE_Data(1.4142_LIBMUSCLE_real8)      ! 64-bit real

! constant list
d8 = LIBMUSCLE_Data_create_list()                   ! empty list
d9 = LIBMUSCLE_Data_create_nils(4_LIBMUSCLE_size)   ! list of length 4
call d9%set_item(1, 'String data')
call d9%set_item(2, .true.)
call d9%set_item(3, 0)
call d9%set_item(4, 1.4142d0)

! dictionary
d10 = LIBMUSCLE_Data_create_dict()
call d10%set_item('x', x)
call d10%set_item('y', y)
call d10%set_item('note', 'Keys must be strings')
call d10%set_item('nest all you want', d9)

! grid
ar = reshape(spread((/1_LIBMUSCLE_int8/), 1, 6), (/2, 3/))
d11 = LIBMUSCLE_Data_create_grid(ar)

! byte array
d12 = LIBMUSCLE_Data_create_byte_array(bytes)

For completeness and as an easy-to-use reference, here’s how to unpack those types:

! check whether a DataConstRef contains a nil value
if (d1%is_nil()) then
    ...
end if

! literals, logicals
character(len=:), allocatable :: c
logical :: l

c = d2%as_character()
l = d3%as_logical()

! various kinds of numbers
integer :: i
integer (int8) :: u8, i8
real (LIBMUSCLE_real4) :: r4
real (LIBMUSCLE_real8) :: r8

i = d4%as_int()
! Fortran has no unsigned numbers, so this is signed
u8 = d5%as_int8()
i8 = d6%as_int8()
r4 = d7%as_real4()
r8 = d8%as_real8()
type(LIBMUSCLE_DataConstRef) :: d9j

! accessing a heterogeneous list
integer (LIBMUSCLE_size) :: j

if (.not. d9%is_a_list()) then
    print *, 'Expected a list'
end if

do j = 1, d9%size()
    d9j = d9%get_item(j)
    if (d9j%is_a_character())
        s = d9j%as_character()
    else if (d9j%is_a_logical())
        l = d9j%as_logical()
    else if (d9j%is_a_int()
        i = d9j%as_int()
    else if (d9j%is_a_real8())
        r8 = d9j%as_real8()
    end if

    call LIBMUSCLE_DataConstRef_free(d9j)
end do

! accessing a dictionary value
type(LIBMUSCLE_DataConstRef) :: d10note
character(len=:), allocatable :: c

if (.not. d10%is_a_dict()) then
    print *, 'Expected a dict'
end if

d10note = d10%get_item('note')
c = d10%as_character()
call LIBMUSCLE_DataConstRef_free(d10note)

! iterating through a dictionary
integer (LIBMUSCLE_size) :: j
character(len=:), allocatable :: k
type(LIBMUSCLE_DataConstRef) :: v

do j = 1, d10%size()
    k = d10%key(j)
    v = d10%value(j)
    ! process here, not shown
    call LIBMUSCLE_DataConstRef_free(v)
end do

! grid
integer (LIBMUSCLE_size), dimension(7) :: shp
integer (LIBMUSCLE_size) :: num_dims
integer (LIBMUSCLE_int8), dimension(2, 3) :: ar
integer :: i
character(len=:), allocatable :: c

if (d11%is_grid_of_real8()) then
    call d11%shape(shp)
    call d11%num_dims(num_dims)
    if (num_dims == 2_LIBMUSCLE_size) then
        call d11%elements(ar)
        c = d11%index(1_LIBMUSCLE_size)
        do i = 1, shp(1)
            print *, c, i, ar(i,:)
        end do
    end if
end if

! byte array
character(len=1), dimension(:), allocatable :: buf
if (d12%is_a_byte_array()) then
    allocate(buf(d12%size()))
    call d12%as_byte_array(buf)
end if

As you can see, sending complex data types with MUSCLE is a bit more difficult in Fortran than in Python, but it is not too burdensome. And Fortran is really good with grids/arrays.

If you want to send large amounts of data, then use grids (arrays) as much as possible. Lists and dicts are more flexible, but they are also slower and take up much more memory. In particular, a set of objects (agents for example) is better sent as a dict of arrays (with one 1D array for each attribute) than as a list of dicts. (The latter will probably work up to 1 million objects or so, but it will still be slower.)

Handling errors

In some cases when you’re asking MUSCLE3 to do something, things can go wrong. For example, receiving a message from a port does not work if the port is not connected (unless you passed a default message, see LIBMUSCLE_Instance_receive()). You cannot extract a real number from a LIBMUSCLE_DataConstRef which contains an integer, or obtain a setting as a character string if the user set it to a logical value. There are two ways of dealing with these situations: checking in advance if it’s going to work and acting accordingly, or just trying and checking whether you were successful. If you do neither (which is easiest), MUSCLE3 will print an error message and stop your program if an error occurs (which is actually usually what you want).

How to check in advance whether an operation can work depends on what you want to do. If you want to receive a message, you could use LIBMUSCLE_Instance_is_connected() to check that the port is connected before trying to receive. For settings that can be of different types, you can use LIBMUSCLE_Instance_is_setting_a_<type>, and for data there is LIBMUSCLE_DataConstRef_is_a_<type>.

This does not work in all cases however, for example if you have received a dictionary that may or may not have a particular key. In that case, (and anywhere else you think it’s more convenient) it’s better to just try, and handle any errors if they occur. This is done by adding one or two extra arguments to the function call that you expect may fail, like this:

type(LIBMUSCLE_DataConstRef) :: rdata, item
integer :: err_code
character(len=:), allocatable :: err_msg

logical :: value

! Receiving rdata omitted

! If there is a key 'key' and its value is a logical, set variable
! value to that, otherwise set it to .true..
item = rdata%get_item('key', err_code, err_msg)
if (err_code /= LIBMUSCLE_success) then
    print *, err_msg
    ! Need to deallocate the message
    deallocate(err_msg)
    value = .true.
else
    value = item%as_logical(err_code)
    if (err_code /= LIBMUSCLE_success) then
        value = .true.
    end if
end if

Note that if an error occurs, an error message will be set and the variable must be deallocated after you’re done using it. If no error occurs, the variable will remain unset, and no deallocation is needed. Passing a variable for an error message is optional, you can also only request the error code, as shown. You cannot get only a message, and no code.

Valid error code values are LIBMUSCLE_success (no error), LIBMUSCLE_domain_error (incorrect input to a function), LIBMUSCLE_out_of_range (invalid index), LIBMUSCLE_logic_error (you tried to do something that doesn’t make sense), LIBMUSCLE_runtime_error (general error, e.g. you’re trying to get the length of a scalar port, which does not make sense), LIBMUSCLE_bad_cast (value doesn’t have the type you’re trying to retrieve it as). See the API Documentation for Fortran for the function you are calling to see which error code it will return when exactly.

(Implementation note: which error is returned when is a bit messy still, and will be cleaned up in a future version.)

Diffusion model

If you’ve studied the above carefully, and have seen the Python version of the diffusion model, then you should now be able to understand the Fortran diffusion model below:

docs/source/examples/fortran/diffusion.f90

! A simple diffusion model on a 1D grid.
!
! The state of this model is a 1D grid of concentrations. It sends out the
! state on each timestep on 'state_out', and can receive an updated state
! on 'state_in' at each state update.

program diffusion
    use ymmsl
    use libmuscle
    implicit none

    type(LIBMUSCLE_PortsDescription) :: ports
    type(LIBMUSCLE_Instance) :: instance

    type(LIBMUSCLE_Message) :: rmsg
    type(LIBMUSCLE_DataConstRef) :: rdata, item


    type(LIBMUSCLE_Message) :: smsg
    type(LIBMUSCLE_Data) :: sdata

    real (selected_real_kind(15)) :: t_cur, t_next, t_max, dt, x_max, dx, d
    integer (LIBMUSCLE_size) :: U_size, n_steps, iteration
    real (selected_real_kind(15)), dimension(:), allocatable :: U, dU
    real (selected_real_kind(15)), dimension(:, :), allocatable :: Us


    ports = LIBMUSCLE_PortsDescription()
    call ports%add(YMMSL_Operator_O_I, 'state_out')
    call ports%add(YMMSL_Operator_S, 'state_in')
    call ports%add(YMMSL_Operator_O_F, 'final_state_out')
    instance = LIBMUSCLE_Instance(ports)
    call LIBMUSCLE_PortsDescription_free(ports)

    do while (instance%reuse_instance())
        ! F_INIT
        t_max = instance%get_setting_as_real8('t_max')
        dt = instance%get_setting_as_real8('dt')
        x_max = instance%get_setting_as_real8('x_max')
        dx = instance%get_setting_as_real8('dx')
        d = instance%get_setting_as_real8('d')

        U_size = nint(x_max / dx)
        allocate (U(U_size), dU(U_size))
        U = 1e-20
        U(26) = 2.0
        U(51) = 2.0
        U(76) = 2.0

        n_steps = int(t_max / dt)
        allocate (Us(U_size, n_steps))

        iteration = 1
        Us(:, iteration) = U

        t_cur = 0.0
        do while (t_cur + dt < t_max .and. iteration < n_steps - 1)
            print *, 't_cur: ', t_cur, 't_max: ', t_max
            ! O_I
            sdata = LIBMUSCLE_Data_create_grid(U, 'x')
            smsg = LIBMUSCLE_Message(t_cur, sdata)
            call LIBMUSCLE_Data_free(sdata)
            t_next = t_cur + dt
            if (t_next + dt <= t_max) then
                call smsg%set_next_timestamp(t_next)
            end if
            call instance%send('state_out', smsg)

            ! S
            rmsg = instance%receive('state_in', smsg)
            rdata = rmsg%get_data()
            call rdata%elements(U)
            call LIBMUSCLE_DataConstRef_free(rdata)
            call LIBMUSCLE_Message_free(rmsg)
            call LIBMUSCLE_Message_free(smsg)

            dU(2:U_size-1) = d * laplacian(U, dx) * dt
            dU(1) = dU(2)
            dU(U_size) = dU(U_size - 1)

            U = U + dU
            iteration = iteration + 1
            Us(:, iteration) = U

            t_cur = t_cur + dt
        end do

        ! O_F
        sdata = LIBMUSCLE_Data_create_grid(U, 'x')
        smsg = LIBMUSCLE_Message(t_cur, sdata)
        call instance%send('final_state_out', smsg)
        call LIBMUSCLE_Message_free(smsg)
        call LIBMUSCLE_Data_free(sdata)
        deallocate (U, dU, Us)
        print *, 'All done'
    end do

    call LIBMUSCLE_Instance_free(instance)

contains

    ! Calculates the Laplacian of vector Z.
    !
    ! @param Z A vector representing a series of samples along a line.
    ! @param dx The spacing between the samples.

    function laplacian(Z, dx)
        real (selected_real_kind(15)), dimension(:), intent(in) :: Z
        real (selected_real_kind(15)), intent(in) :: dx
        real (selected_real_kind(15)), allocatable, dimension(:) :: laplacian
        integer :: n

        n = size(Z)
        allocate(laplacian(size(Z) - 2))
        laplacian = (Z(1:n-2) + Z(3:n) - 2.0d0 * Z(2:n-1)) / (dx * dx)

    end function laplacian

end program diffusion

In docs/source/examples/fortran you will find Fortran implementations of the reaction and diffusion models, as well as of the Monte Carlo sampler and round-robin load balancer. See the Uncertainty Quantification section to learn what those are for.