Repository: multiscale/muscle3 (language: python)

Issue:
bool8 deprecated in numpy
In NumPy 1.24, `numpy.bool8` is deprecated, and just to make sure no one misses it, it's already been removed from the type annotations, so now mypy fails. The new `bool8` is `bool_`, which has been around since at least 1.13, so we can safely replace it.

I guess someone hit the 8 with their software carpentry hammer?


Discussion / hints:

Fix patch:
diff --git a/Doxyfile b/Doxyfile
index f0c512dd..a8b7d7b2 100644
--- a/Doxyfile
+++ b/Doxyfile
@@ -1,4 +1,4 @@
-# Doxyfile 1.8.13
+# Doxyfile 1.8.17
 
 # This file describes the settings to be used by the documentation system
 # doxygen (www.doxygen.org) for a project.
@@ -17,11 +17,11 @@
 # Project related configuration options
 #---------------------------------------------------------------------------
 
-# This tag specifies the encoding used for all characters in the config file
-# that follow. The default is UTF-8 which is also the encoding used for all text
-# before the first occurrence of this tag. Doxygen uses libiconv (or the iconv
-# built into libc) for the transcoding. See http://www.gnu.org/software/libiconv
-# for the list of possible encodings.
+# This tag specifies the encoding used for all characters in the configuration
+# file that follow. The default is UTF-8 which is also the encoding used for all
+# text before the first occurrence of this tag. Doxygen uses libiconv (or the
+# iconv built into libc) for the transcoding. See
+# https://www.gnu.org/software/libiconv/ for the list of possible encodings.
 # The default value is: UTF-8.
 
 DOXYFILE_ENCODING      = UTF-8
@@ -93,6 +93,14 @@ ALLOW_UNICODE_NAMES    = NO
 
 OUTPUT_LANGUAGE        = English
 
+# The OUTPUT_TEXT_DIRECTION tag is used to specify the direction in which all
+# documentation generated by doxygen is written. Doxygen will use this
+# information to generate all generated output in the proper direction.
+# Possible values are: None, LTR, RTL and Context.
+# The default value is: None.
+
+OUTPUT_TEXT_DIRECTION  = None
+
 # If the BRIEF_MEMBER_DESC tag is set to YES, doxygen will include brief member
 # descriptions after the members that are listed in the file and class
 # documentation (similar to Javadoc). Set to NO to disable this.
@@ -189,6 +197,16 @@ SHORT_NAMES            = NO
 
 JAVADOC_AUTOBRIEF      = YES
 
+# If the JAVADOC_BANNER tag is set to YES then doxygen will interpret a line
+# such as
+# /***************
+# as being the beginning of a Javadoc-style comment "banner". If set to NO, the
+# Javadoc-style will behave just like regular comments and it will not be
+# interpreted by doxygen.
+# The default value is: NO.
+
+JAVADOC_BANNER         = NO
+
 # If the QT_AUTOBRIEF tag is set to YES then doxygen will interpret the first
 # line (until the first dot) of a Qt-style comment as the brief description. If
 # set to NO, the Qt-style will behave just like regular Qt-style comments (thus
@@ -236,7 +254,12 @@ TAB_SIZE               = 4
 # will allow you to put the command \sideeffect (or @sideeffect) in the
 # documentation, which will result in a user-defined paragraph with heading
 # "Side Effects:". You can put \n's in the value part of an alias to insert
-# newlines.
+# newlines (in the resulting output). You can put ^^ in the value part of an
+# alias to insert a newline as if a physical newline was in the original file.
+# When you need a literal { or } or , in the value part of an alias you have to
+# escape them by means of a backslash (\), this can lead to conflicts with the
+# commands \{ and \} for these it is advised to use the version @{ and @} or use
+# a double escape (\\{ and \\})
 
 ALIASES                =
 
@@ -274,17 +297,26 @@ OPTIMIZE_FOR_FORTRAN   = NO
 
 OPTIMIZE_OUTPUT_VHDL   = NO
 
+# Set the OPTIMIZE_OUTPUT_SLICE tag to YES if your project consists of Slice
+# sources only. Doxygen will then generate output that is more tailored for that
+# language. For instance, namespaces will be presented as modules, types will be
+# separated into more groups, etc.
+# The default value is: NO.
+
+OPTIMIZE_OUTPUT_SLICE  = NO
+
 # Doxygen selects the parser to use depending on the extension of the files it
 # parses. With this tag you can assign which parser to use for a given
 # extension. Doxygen has a built-in mapping, but you can override or extend it
 # using this tag. The format is ext=language, where ext is a file extension, and
-# language is one of the parsers supported by doxygen: IDL, Java, Javascript,
-# C#, C, C++, D, PHP, Objective-C, Python, Fortran (fixed format Fortran:
-# FortranFixed, free formatted Fortran: FortranFree, unknown formatted Fortran:
-# Fortran. In the later case the parser tries to guess whether the code is fixed
-# or free formatted code, this is the default for Fortran type files), VHDL. For
-# instance to make doxygen treat .inc files as Fortran files (default is PHP),
-# and .f files as C (default is Fortran), use: inc=Fortran f=C.
+# language is one of the parsers supported by doxygen: IDL, Java, JavaScript,
+# Csharp (C#), C, C++, D, PHP, md (Markdown), Objective-C, Python, Slice,
+# Fortran (fixed format Fortran: FortranFixed, free formatted Fortran:
+# FortranFree, unknown formatted Fortran: Fortran. In the later case the parser
+# tries to guess whether the code is fixed or free formatted code, this is the
+# default for Fortran type files), VHDL, tcl. For instance to make doxygen treat
+# .inc files as Fortran files (default is PHP), and .f files as C (default is
+# Fortran), use: inc=Fortran f=C.
 #
 # Note: For files without extension you can use no_extension as a placeholder.
 #
@@ -295,7 +327,7 @@ EXTENSION_MAPPING      =
 
 # If the MARKDOWN_SUPPORT tag is enabled then doxygen pre-processes all comments
 # according to the Markdown format, which allows for more readable
-# documentation. See http://daringfireball.net/projects/markdown/ for details.
+# documentation. See https://daringfireball.net/projects/markdown/ for details.
 # The output of markdown processing is further processed by doxygen, so you can
 # mix doxygen, HTML, and XML commands with Markdown formatting. Disable only in
 # case of backward compatibilities issues.
@@ -307,7 +339,7 @@ MARKDOWN_SUPPORT       = YES
 # to that level are automatically included in the table of contents, even if
 # they do not have an id attribute.
 # Note: This feature currently applies only to Markdown headings.
-# Minimum value: 0, maximum value: 99, default value: 0.
+# Minimum value: 0, maximum value: 99, default value: 5.
 # This tag requires that the tag MARKDOWN_SUPPORT is set to YES.
 
 TOC_INCLUDE_HEADINGS   = 0
@@ -337,7 +369,7 @@ BUILTIN_STL_SUPPORT    = NO
 CPP_CLI_SUPPORT        = NO
 
 # Set the SIP_SUPPORT tag to YES if your project consists of sip (see:
-# http://www.riverbankcomputing.co.uk/software/sip/intro) sources only. Doxygen
+# https://www.riverbankcomputing.com/software/sip/intro) sources only. Doxygen
 # will parse them like normal C++ but will assume all classes use public instead
 # of private inheritance when no explicit protection keyword is present.
 # The default value is: NO.
@@ -443,6 +475,12 @@ EXTRACT_ALL            = NO
 
 EXTRACT_PRIVATE        = NO
 
+# If the EXTRACT_PRIV_VIRTUAL tag is set to YES, documented private virtual
+# methods of a class will be included in the documentation.
+# The default value is: NO.
+
+EXTRACT_PRIV_VIRTUAL   = NO
+
 # If the EXTRACT_PACKAGE tag is set to YES, all members with package or internal
 # scope will be included in the documentation.
 # The default value is: NO.
@@ -497,8 +535,8 @@ HIDE_UNDOC_MEMBERS     = NO
 HIDE_UNDOC_CLASSES     = NO
 
 # If the HIDE_FRIEND_COMPOUNDS tag is set to YES, doxygen will hide all friend
-# (class|struct|union) declarations. If set to NO, these declarations will be
-# included in the documentation.
+# declarations. If set to NO, these declarations will be included in the
+# documentation.
 # The default value is: NO.
 
 HIDE_FRIEND_COMPOUNDS  = YES
@@ -521,7 +559,7 @@ INTERNAL_DOCS          = NO
 # names in lower-case letters. If set to YES, upper-case letters are also
 # allowed. This is useful if you have classes or files whose names only differ
 # in case and if your file system supports case sensitive file names. Windows
-# and Mac users are advised to set this option to NO.
+# (including Cygwin) ands Mac users are advised to set this option to NO.
 # The default value is: system dependent.
 
 CASE_SENSE_NAMES       = YES
@@ -708,7 +746,7 @@ LAYOUT_FILE            =
 # The CITE_BIB_FILES tag can be used to specify one or more bib files containing
 # the reference definitions. This must be a list of .bib files. The .bib
 # extension is automatically appended if omitted. This requires the bibtex tool
-# to be installed. See also http://en.wikipedia.org/wiki/BibTeX for more info.
+# to be installed. See also https://en.wikipedia.org/wiki/BibTeX for more info.
 # For LaTeX the style of the bibliography can be controlled using
 # LATEX_BIB_STYLE. To use this feature you need bibtex and perl available in the
 # search path. See also \cite for info how to create references.
@@ -753,7 +791,8 @@ WARN_IF_DOC_ERROR      = YES
 # This WARN_NO_PARAMDOC option can be enabled to get warnings for functions that
 # are documented, but have no documentation for their parameters or return
 # value. If set to NO, doxygen will only warn about wrong or incomplete
-# parameter documentation, but not about the absence of documentation.
+# parameter documentation, but not about the absence of documentation. If
+# EXTRACT_ALL is set to YES then this flag will automatically be disabled.
 # The default value is: NO.
 
 WARN_NO_PARAMDOC       = NO
@@ -795,7 +834,7 @@ INPUT                  = libmuscle/cpp/src
 # This tag can be used to specify the character encoding of the source files
 # that doxygen parses. Internally doxygen uses the UTF-8 encoding. Doxygen uses
 # libiconv (or the iconv built into libc) for the transcoding. See the libiconv
-# documentation (see: http://www.gnu.org/software/libiconv) for the list of
+# documentation (see: https://www.gnu.org/software/libiconv/) for the list of
 # possible encodings.
 # The default value is: UTF-8.
 
@@ -812,8 +851,10 @@ INPUT_ENCODING         = UTF-8
 # If left blank the following patterns are tested:*.c, *.cc, *.cxx, *.cpp,
 # *.c++, *.java, *.ii, *.ixx, *.ipp, *.i++, *.inl, *.idl, *.ddl, *.odl, *.h,
 # *.hh, *.hxx, *.hpp, *.h++, *.cs, *.d, *.php, *.php4, *.php5, *.phtml, *.inc,
-# *.m, *.markdown, *.md, *.mm, *.dox, *.py, *.pyw, *.f90, *.f95, *.f03, *.f08,
-# *.f, *.for, *.tcl, *.vhd, *.vhdl, *.ucf and *.qsf.
+# *.m, *.markdown, *.md, *.mm, *.dox (to be provided as doxygen C comment),
+# *.doc (to be provided as doxygen C comment), *.txt (to be provided as doxygen
+# C comment), *.py, *.pyw, *.f90, *.f95, *.f03, *.f08, *.f, *.for, *.tcl, *.vhd,
+# *.vhdl, *.ucf, *.qsf and *.ice.
 
 FILE_PATTERNS          = *.c \
                          *.cc \
@@ -873,7 +914,10 @@ RECURSIVE              = YES
 # Note that relative paths are relative to the directory from which doxygen is
 # run.
 
-EXCLUDE                = libmuscle/cpp/src/muscle_manager_protocol
+EXCLUDE                = libmuscle/cpp/src/libmuscle/bindings \
+                         libmuscle/cpp/src/libmuscle/tests \
+                         libmuscle/cpp/src/ymmsl/bindings \
+                         libmuscle/cpp/src/ymmsl/tests
 
 # The EXCLUDE_SYMLINKS tag can be used to select whether or not files or
 # directories that are symbolic links (a Unix file system feature) are excluded
@@ -1011,7 +1055,7 @@ INLINE_SOURCES         = NO
 STRIP_CODE_COMMENTS    = YES
 
 # If the REFERENCED_BY_RELATION tag is set to YES then for each documented
-# function all documented functions referencing it will be listed.
+# entity all documented functions referencing it will be listed.
 # The default value is: NO.
 
 REFERENCED_BY_RELATION = NO
@@ -1043,12 +1087,12 @@ SOURCE_TOOLTIPS        = YES
 # If the USE_HTAGS tag is set to YES then the references to source code will
 # point to the HTML generated by the htags(1) tool instead of doxygen built-in
 # source browser. The htags tool is part of GNU's global source tagging system
-# (see http://www.gnu.org/software/global/global.html). You will need version
+# (see https://www.gnu.org/software/global/global.html). You will need version
 # 4.8.6 or higher.
 #
 # To use it do the following:
 # - Install the latest version of global
-# - Enable SOURCE_BROWSER and USE_HTAGS in the config file
+# - Enable SOURCE_BROWSER and USE_HTAGS in the configuration file
 # - Make sure the INPUT points to the root of the source tree
 # - Run doxygen as normal
 #
@@ -1076,7 +1120,7 @@ VERBATIM_HEADERS       = YES
 # rich C++ code for which doxygen's built-in parser lacks the necessary type
 # information.
 # Note: The availability of this option depends on whether or not doxygen was
-# generated with the -Duse-libclang=ON option for CMake.
+# generated with the -Duse_libclang=ON option for CMake.
 # The default value is: NO.
 
 CLANG_ASSISTED_PARSING = NO
@@ -1089,6 +1133,16 @@ CLANG_ASSISTED_PARSING = NO
 
 CLANG_OPTIONS          =
 
+# If clang assisted parsing is enabled you can provide the clang parser with the
+# path to the compilation database (see:
+# http://clang.llvm.org/docs/HowToSetupToolingForLLVM.html) used when the files
+# were built. This is equivalent to specifying the "-p" option to a clang tool,
+# such as clang-check. These options will then be passed to the parser.
+# Note: The availability of this option depends on whether or not doxygen was
+# generated with the -Duse_libclang=ON option for CMake.
+
+CLANG_DATABASE_PATH    =
+
 #---------------------------------------------------------------------------
 # Configuration options related to the alphabetical class index
 #---------------------------------------------------------------------------
@@ -1122,7 +1176,7 @@ IGNORE_PREFIX          =
 # If the GENERATE_HTML tag is set to YES, doxygen will generate HTML output
 # The default value is: YES.
 
-GENERATE_HTML          = YES
+GENERATE_HTML          = NO
 
 # The HTML_OUTPUT tag is used to specify where the HTML docs will be put. If a
 # relative path is entered the value of OUTPUT_DIRECTORY will be put in front of
@@ -1207,7 +1261,7 @@ HTML_EXTRA_FILES       =
 # The HTML_COLORSTYLE_HUE tag controls the color of the HTML output. Doxygen
 # will adjust the colors in the style sheet and background images according to
 # this color. Hue is specified as an angle on a colorwheel, see
-# http://en.wikipedia.org/wiki/Hue for more information. For instance the value
+# https://en.wikipedia.org/wiki/Hue for more information. For instance the value
 # 0 represents red, 60 is yellow, 120 is green, 180 is cyan, 240 is blue, 300
 # purple, and 360 is red again.
 # Minimum value: 0, maximum value: 359, default value: 220.
@@ -1243,6 +1297,17 @@ HTML_COLORSTYLE_GAMMA  = 80
 
 HTML_TIMESTAMP         = NO
 
+# If the HTML_DYNAMIC_MENUS tag is set to YES then the generated HTML
+# documentation will contain a main index with vertical navigation menus that
+# are dynamically created via JavaScript. If disabled, the navigation index will
+# consists of multiple levels of tabs that are statically embedded in every HTML
+# page. Disable this option to support browsers that do not have JavaScript,
+# like the Qt help browser.
+# The default value is: YES.
+# This tag requires that the tag GENERATE_HTML is set to YES.
+
+HTML_DYNAMIC_MENUS     = YES
+
 # If the HTML_DYNAMIC_SECTIONS tag is set to YES then the generated HTML
 # documentation will contain sections that can be hidden and shown after the
 # page has loaded.
@@ -1266,13 +1331,13 @@ HTML_INDEX_NUM_ENTRIES = 100
 
 # If the GENERATE_DOCSET tag is set to YES, additional index files will be
 # generated that can be used as input for Apple's Xcode 3 integrated development
-# environment (see: http://developer.apple.com/tools/xcode/), introduced with
-# OSX 10.5 (Leopard). To create a documentation set, doxygen will generate a
+# environment (see: https://developer.apple.com/xcode/), introduced with OSX
+# 10.5 (Leopard). To create a documentation set, doxygen will generate a
 # Makefile in the HTML output directory. Running make will produce the docset in
 # that directory and running make install will install the docset in
 # ~/Library/Developer/Shared/Documentation/DocSets so that Xcode will find it at
-# startup. See http://developer.apple.com/tools/creatingdocsetswithdoxygen.html
-# for more information.
+# startup. See https://developer.apple.com/library/archive/featuredarticles/Doxy
+# genXcode/_index.html for more information.
 # The default value is: NO.
 # This tag requires that the tag GENERATE_HTML is set to YES.
 
@@ -1311,7 +1376,7 @@ DOCSET_PUBLISHER_NAME  = Publisher
 # If the GENERATE_HTMLHELP tag is set to YES then doxygen generates three
 # additional HTML index files: index.hhp, index.hhc, and index.hhk. The
 # index.hhp is a project file that can be read by Microsoft's HTML Help Workshop
-# (see: http://www.microsoft.com/en-us/download/details.aspx?id=21138) on
+# (see: https://www.microsoft.com/en-us/download/details.aspx?id=21138) on
 # Windows.
 #
 # The HTML Help Workshop contains a compiler that can convert all HTML output
@@ -1387,7 +1452,7 @@ QCH_FILE               =
 
 # The QHP_NAMESPACE tag specifies the namespace to use when generating Qt Help
 # Project output. For more information please see Qt Help Project / Namespace
-# (see: http://qt-project.org/doc/qt-4.8/qthelpproject.html#namespace).
+# (see: https://doc.qt.io/archives/qt-4.8/qthelpproject.html#namespace).
 # The default value is: org.doxygen.Project.
 # This tag requires that the tag GENERATE_QHP is set to YES.
 
@@ -1395,7 +1460,7 @@ QHP_NAMESPACE          = org.doxygen.Project
 
 # The QHP_VIRTUAL_FOLDER tag specifies the namespace to use when generating Qt
 # Help Project output. For more information please see Qt Help Project / Virtual
-# Folders (see: http://qt-project.org/doc/qt-4.8/qthelpproject.html#virtual-
+# Folders (see: https://doc.qt.io/archives/qt-4.8/qthelpproject.html#virtual-
 # folders).
 # The default value is: doc.
 # This tag requires that the tag GENERATE_QHP is set to YES.
@@ -1404,7 +1469,7 @@ QHP_VIRTUAL_FOLDER     = doc
 
 # If the QHP_CUST_FILTER_NAME tag is set, it specifies the name of a custom
 # filter to add. For more information please see Qt Help Project / Custom
-# Filters (see: http://qt-project.org/doc/qt-4.8/qthelpproject.html#custom-
+# Filters (see: https://doc.qt.io/archives/qt-4.8/qthelpproject.html#custom-
 # filters).
 # This tag requires that the tag GENERATE_QHP is set to YES.
 
@@ -1412,7 +1477,7 @@ QHP_CUST_FILTER_NAME   =
 
 # The QHP_CUST_FILTER_ATTRS tag specifies the list of the attributes of the
 # custom filter to add. For more information please see Qt Help Project / Custom
-# Filters (see: http://qt-project.org/doc/qt-4.8/qthelpproject.html#custom-
+# Filters (see: https://doc.qt.io/archives/qt-4.8/qthelpproject.html#custom-
 # filters).
 # This tag requires that the tag GENERATE_QHP is set to YES.
 
@@ -1420,7 +1485,7 @@ QHP_CUST_FILTER_ATTRS  =
 
 # The QHP_SECT_FILTER_ATTRS tag specifies the list of the attributes this
 # project's filter section matches. Qt Help Project / Filter Attributes (see:
-# http://qt-project.org/doc/qt-4.8/qthelpproject.html#filter-attributes).
+# https://doc.qt.io/archives/qt-4.8/qthelpproject.html#filter-attributes).
 # This tag requires that the tag GENERATE_QHP is set to YES.
 
 QHP_SECT_FILTER_ATTRS  =
@@ -1513,7 +1578,7 @@ EXT_LINKS_IN_WINDOW    = NO
 
 FORMULA_FONTSIZE       = 10
 
-# Use the FORMULA_TRANPARENT tag to determine whether or not the images
+# Use the FORMULA_TRANSPARENT tag to determine whether or not the images
 # generated for formulas are transparent PNGs. Transparent PNGs are not
 # supported properly for IE 6.0, but are supported on all modern browsers.
 #
@@ -1524,8 +1589,14 @@ FORMULA_FONTSIZE       = 10
 
 FORMULA_TRANSPARENT    = YES
 
+# The FORMULA_MACROFILE can contain LaTeX \newcommand and \renewcommand commands
+# to create new LaTeX commands to be used in formulas as building blocks. See
+# the section "Including formulas" for details.
+
+FORMULA_MACROFILE      =
+
 # Enable the USE_MATHJAX option to render LaTeX formulas using MathJax (see
-# http://www.mathjax.org) which uses client side Javascript for the rendering
+# https://www.mathjax.org) which uses client side JavaScript for the rendering
 # instead of using pre-rendered bitmaps. Use this if you do not have LaTeX
 # installed or if you want to formulas look prettier in the HTML output. When
 # enabled you may also need to install MathJax separately and configure the path
@@ -1552,8 +1623,8 @@ MATHJAX_FORMAT         = HTML-CSS
 # MATHJAX_RELPATH should be ../mathjax. The default value points to the MathJax
 # Content Delivery Network so you can quickly see the result without installing
 # MathJax. However, it is strongly recommended to install a local copy of
-# MathJax from http://www.mathjax.org before deployment.
-# The default value is: http://cdn.mathjax.org/mathjax/latest.
+# MathJax from https://www.mathjax.org before deployment.
+# The default value is: https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.5/.
 # This tag requires that the tag USE_MATHJAX is set to YES.
 
 MATHJAX_RELPATH        = http://cdn.mathjax.org/mathjax/latest
@@ -1595,7 +1666,7 @@ MATHJAX_CODEFILE       =
 SEARCHENGINE           = YES
 
 # When the SERVER_BASED_SEARCH tag is enabled the search engine will be
-# implemented using a web server instead of a web client using Javascript. There
+# implemented using a web server instead of a web client using JavaScript. There
 # are two flavors of web server based searching depending on the EXTERNAL_SEARCH
 # setting. When disabled, doxygen will generate a PHP script for searching and
 # an index file used by the script. When EXTERNAL_SEARCH is enabled the indexing
@@ -1614,7 +1685,7 @@ SERVER_BASED_SEARCH    = NO
 #
 # Doxygen ships with an example indexer (doxyindexer) and search engine
 # (doxysearch.cgi) which are based on the open source search engine library
-# Xapian (see: http://xapian.org/).
+# Xapian (see: https://xapian.org/).
 #
 # See the section "External Indexing and Searching" for details.
 # The default value is: NO.
@@ -1627,7 +1698,7 @@ EXTERNAL_SEARCH        = NO
 #
 # Doxygen ships with an example indexer (doxyindexer) and search engine
 # (doxysearch.cgi) which are based on the open source search engine library
-# Xapian (see: http://xapian.org/). See the section "External Indexing and
+# Xapian (see: https://xapian.org/). See the section "External Indexing and
 # Searching" for details.
 # This tag requires that the tag SEARCHENGINE is set to YES.
 
@@ -1679,21 +1750,35 @@ LATEX_OUTPUT           = latex
 # The LATEX_CMD_NAME tag can be used to specify the LaTeX command name to be
 # invoked.
 #
-# Note that when enabling USE_PDFLATEX this option is only used for generating
-# bitmaps for formulas in the HTML output, but not in the Makefile that is
-# written to the output directory.
-# The default file is: latex.
+# Note that when not enabling USE_PDFLATEX the default is latex when enabling
+# USE_PDFLATEX the default is pdflatex and when in the later case latex is
+# chosen this is overwritten by pdflatex. For specific output languages the
+# default can have been set differently, this depends on the implementation of
+# the output language.
 # This tag requires that the tag GENERATE_LATEX is set to YES.
 
 LATEX_CMD_NAME         = latex
 
 # The MAKEINDEX_CMD_NAME tag can be used to specify the command name to generate
 # index for LaTeX.
+# Note: This tag is used in the Makefile / make.bat.
+# See also: LATEX_MAKEINDEX_CMD for the part in the generated output file
+# (.tex).
 # The default file is: makeindex.
 # This tag requires that the tag GENERATE_LATEX is set to YES.
 
 MAKEINDEX_CMD_NAME     = makeindex
 
+# The LATEX_MAKEINDEX_CMD tag can be used to specify the command name to
+# generate index for LaTeX. In case there is no backslash (\) as first character
+# it will be automatically added in the LaTeX code.
+# Note: This tag is used in the generated output file (.tex).
+# See also: MAKEINDEX_CMD_NAME for the part in the Makefile / make.bat.
+# The default value is: makeindex.
+# This tag requires that the tag GENERATE_LATEX is set to YES.
+
+LATEX_MAKEINDEX_CMD    = makeindex
+
 # If the COMPACT_LATEX tag is set to YES, doxygen generates more compact LaTeX
 # documents. This may be useful for small projects and may help to save some
 # trees in general.
@@ -1814,7 +1899,7 @@ LATEX_SOURCE_CODE      = NO
 
 # The LATEX_BIB_STYLE tag can be used to specify the style to use for the
 # bibliography, e.g. plainnat, or ieeetr. See
-# http://en.wikipedia.org/wiki/BibTeX and \cite for more info.
+# https://en.wikipedia.org/wiki/BibTeX and \cite for more info.
 # The default value is: plain.
 # This tag requires that the tag GENERATE_LATEX is set to YES.
 
@@ -1828,6 +1913,14 @@ LATEX_BIB_STYLE        = plain
 
 LATEX_TIMESTAMP        = NO
 
+# The LATEX_EMOJI_DIRECTORY tag is used to specify the (relative or absolute)
+# path from which the emoji images will be read. If a relative path is entered,
+# it will be relative to the LATEX_OUTPUT directory. If left blank the
+# LATEX_OUTPUT directory will be used.
+# This tag requires that the tag GENERATE_LATEX is set to YES.
+
+LATEX_EMOJI_DIRECTORY  =
+
 #---------------------------------------------------------------------------
 # Configuration options related to the RTF output
 #---------------------------------------------------------------------------
@@ -1867,9 +1960,9 @@ COMPACT_RTF            = NO
 
 RTF_HYPERLINKS         = NO
 
-# Load stylesheet definitions from file. Syntax is similar to doxygen's config
-# file, i.e. a series of assignments. You only have to provide replacements,
-# missing definitions are set to their default value.
+# Load stylesheet definitions from file. Syntax is similar to doxygen's
+# configuration file, i.e. a series of assignments. You only have to provide
+# replacements, missing definitions are set to their default value.
 #
 # See also section "Doxygen usage" for information on how to generate the
 # default style sheet that doxygen normally uses.
@@ -1878,8 +1971,8 @@ RTF_HYPERLINKS         = NO
 RTF_STYLESHEET_FILE    =
 
 # Set optional variables used in the generation of an RTF document. Syntax is
-# similar to doxygen's config file. A template extensions file can be generated
-# using doxygen -e rtf extensionFile.
+# similar to doxygen's configuration file. A template extensions file can be
+# generated using doxygen -e rtf extensionFile.
 # This tag requires that the tag GENERATE_RTF is set to YES.
 
 RTF_EXTENSIONS_FILE    =
@@ -1965,6 +2058,13 @@ XML_OUTPUT             = xml
 
 XML_PROGRAMLISTING     = YES
 
+# If the XML_NS_MEMB_FILE_SCOPE tag is set to YES, doxygen will include
+# namespace members in file scope as well, matching the HTML output.
+# The default value is: NO.
+# This tag requires that the tag GENERATE_XML is set to YES.
+
+XML_NS_MEMB_FILE_SCOPE = NO
+
 #---------------------------------------------------------------------------
 # Configuration options related to the DOCBOOK output
 #---------------------------------------------------------------------------
@@ -1997,9 +2097,9 @@ DOCBOOK_PROGRAMLISTING = NO
 #---------------------------------------------------------------------------
 
 # If the GENERATE_AUTOGEN_DEF tag is set to YES, doxygen will generate an
-# AutoGen Definitions (see http://autogen.sf.net) file that captures the
-# structure of the code including all documentation. Note that this feature is
-# still experimental and incomplete at the moment.
+# AutoGen Definitions (see http://autogen.sourceforge.net/) file that captures
+# the structure of the code including all documentation. Note that this feature
+# is still experimental and incomplete at the moment.
 # The default value is: NO.
 
 GENERATE_AUTOGEN_DEF   = NO
@@ -2099,7 +2199,8 @@ INCLUDE_FILE_PATTERNS  =
 # recursively expanded use the := operator instead of the = operator.
 # This tag requires that the tag ENABLE_PREPROCESSING is set to YES.
 
-PREDEFINED             = MUSCLE_ENABLE_MPI DOXYGEN_SHOULD_SKIP_THIS
+PREDEFINED             = MUSCLE_ENABLE_MPI \
+                         DOXYGEN_SHOULD_SKIP_THIS
 
 # If the MACRO_EXPANSION and EXPAND_ONLY_PREDEF tags are set to YES then this
 # tag can be used to specify a list of macro names that should be expanded. The
@@ -2166,12 +2267,6 @@ EXTERNAL_GROUPS        = YES
 
 EXTERNAL_PAGES         = YES
 
-# The PERL_PATH should be the absolute path and name of the perl script
-# interpreter (i.e. the result of 'which perl').
-# The default file (with absolute path) is: /usr/bin/perl.
-
-PERL_PATH              = /usr/bin/perl
-
 #---------------------------------------------------------------------------
 # Configuration options related to the dot tool
 #---------------------------------------------------------------------------
@@ -2185,15 +2280,6 @@ PERL_PATH              = /usr/bin/perl
 
 CLASS_DIAGRAMS         = YES
 
-# You can define message sequence charts within doxygen comments using the \msc
-# command. Doxygen will then run the mscgen tool (see:
-# http://www.mcternan.me.uk/mscgen/)) to produce the chart and insert it in the
-# documentation. The MSCGEN_PATH tag allows you to specify the directory where
-# the mscgen tool resides. If left empty the tool is assumed to be found in the
-# default search path.
-
-MSCGEN_PATH            =
-
 # You can include diagrams made with dia in doxygen documentation. Doxygen will
 # then run dia to produce the diagram and insert it in the documentation. The
 # DIA_PATH tag allows you to specify the directory where the dia binary resides.
diff --git a/docs/source/conf.py b/docs/source/conf.py
index 7464d1c0..5cef345e 100644
--- a/docs/source/conf.py
+++ b/docs/source/conf.py
@@ -47,6 +47,7 @@
         'breathe',
         'sphinx.ext.autodoc',
         'sphinx.ext.autosectionlabel',
+        'sphinx.ext.intersphinx',
         'sphinx.ext.napoleon',
         'sphinx.ext.todo',
         'sphinx.ext.viewcode',
@@ -109,6 +110,15 @@
 
 breathe_default_members = ('members',)
 
+# Configuration of sphinx.ext.intersphinx
+# See https://www.sphinx-doc.org/en/master/usage/extensions/intersphinx.html
+intersphinx_mapping = {
+    "python": ("https://docs.python.org/3/", None),
+    "numpy": ("https://numpy.org/doc/stable", None),
+    "ymmsl": ("https://ymmsl-python.readthedocs.io/en/stable", None),
+}
+
+
 # -- Patch version into installation instructions --
 def patch_installation_version():
     with open('installing.rst', 'w') as out_file:
diff --git a/docs/source/examples/python/interact_coupling.py b/docs/source/examples/python/interact_coupling.py
index ff9408c7..3df5e11e 100644
--- a/docs/source/examples/python/interact_coupling.py
+++ b/docs/source/examples/python/interact_coupling.py
@@ -1,5 +1,5 @@
 import logging
-from typing import Any, Optional, Tuple
+from typing import Any, Optional, Tuple, Dict
 
 from libmuscle import Instance, Message
 from libmuscle.runner import run_simulation
@@ -129,7 +129,8 @@ class Peer:
     via the instance object.
     """
     def __init__(
-            self, instance: Instance, in_port: str, out_port: str) -> None:
+            self, instance: Instance, in_port: str, out_port: str,
+            resume_from_state: Any = None) -> None:
         """Create a Peer object.
 
         This also receives an initial message from the peer model, and
@@ -145,11 +146,20 @@ def __init__(
         self.out_port = out_port
         self.cache = DataCache()
 
-        msg = self.instance.receive(self.in_port)
-        self.cache.add_data(msg.timestamp, msg.data)
-        self.rcvd = msg.timestamp
-        self.to_send = msg.timestamp
-        self.next = msg.next_timestamp
+        if resume_from_state:
+            self.cache.t_cur = resume_from_state['cache.t_cur']
+            self.cache.data_cur = resume_from_state['cache.data_cur']
+            self.cache.t_next = resume_from_state['cache.t_next']
+            self.cache.data_next = resume_from_state['cache.data_next']
+            self.rcvd = resume_from_state['rcvd']
+            self.to_send = resume_from_state['to_send']
+            self.next = resume_from_state['next']
+        else:
+            msg = self.instance.receive(self.in_port)
+            self.cache.add_data(msg.timestamp, msg.data)
+            self.rcvd = msg.timestamp
+            self.to_send = msg.timestamp
+            self.next = msg.next_timestamp
 
     def done(self) -> bool:
         """Return whether we are done commmunicating with this peer."""
@@ -200,6 +210,17 @@ def send(self, t: float, data: Any) -> None:
         self.instance.send(self.out_port, Message(t, self.next, data))
         self.to_send = self.next
 
+    def get_state(self) -> Dict[str, Any]:
+        """Return the current state of this object as a MUSCLE-serializable dict
+        """
+        return {'cache.t_cur': self.cache.t_cur,
+                'cache.data_cur': self.cache.data_cur,
+                'cache.t_next': self.cache.t_next,
+                'cache.data_next': self.cache.data_next,
+                'rcvd': self.rcvd,
+                'to_send': self.to_send,
+                'next': self.next}
+
 
 def temporal_coupler() -> None:
     """Model component connecting two scale-overlapping submodels.
@@ -241,6 +262,56 @@ def temporal_coupler() -> None:
                 b.send(t, data)
 
 
+def checkpointing_temporal_coupler() -> None:
+    """Model component connecting two scale-overlapping submodels.
+
+    This component sits in between two scale-overlapping submodels
+    running at different (and potentially variable) timesteps and
+    ensures that each of these peers receives a message whenever it
+    expects one, and can send a message whenever it expects to do so.
+
+    This function extends :func:`temporal_coupler` with checkpointing
+    capabilities.
+    """
+    instance = Instance({
+        Operator.O_I: ['a_out', 'b_out'],
+        Operator.S: ['a_in', 'b_in']})
+
+    while instance.reuse_instance():
+        if instance.resuming():
+            state = instance.load_snapshot().data
+            if state is not None:
+                a = Peer(instance, 'a_in', 'a_out', state['a'])
+                b = Peer(instance, 'b_in', 'b_out', state['b'])
+
+        if instance.should_init():
+            # Receive initial messages and initialise state
+            a = Peer(instance, 'a_in', 'a_out')
+            b = Peer(instance, 'b_in', 'b_out')
+
+        # Send and receive as needed
+        while not a.done() or not b.done():
+            if a.can_receive():
+                a.receive()
+            elif b.can_receive():
+                b.receive()
+            elif a.can_send(b.rcvd, b.next):
+                t, data = b.cache.get_data(a.to_send)
+                a.send(t, data)
+            elif b.can_send(a.rcvd, a.next):
+                t, data = a.cache.get_data(b.to_send)
+                b.send(t, data)
+
+            t_cur = min(a.rcvd, b.rcvd)
+            if instance.should_save_snapshot(t_cur):
+                instance.save_snapshot(Message(
+                        t_cur, None, {'a': a.get_state(), 'b': b.get_state()}))
+
+        t_cur = min(a.rcvd, b.rcvd)
+        if instance.should_save_final_snapshot():
+            instance.save_final_snapshot(Message(t_cur, None, None))
+
+
 if __name__ == '__main__':
     logging.basicConfig()
     logging.getLogger().setLevel(logging.INFO)
diff --git a/docs/source/examples/python/requirements.txt b/docs/source/examples/python/requirements.txt
index fa14df52..cabe0c71 100644
--- a/docs/source/examples/python/requirements.txt
+++ b/docs/source/examples/python/requirements.txt
@@ -3,6 +3,6 @@ numpy<1.22; python_version=='3.7'
 numpy>=1.22,<=1.25; python_version>='3.8'
 sobol_seq==0.2.0
 yatiml==0.9.0
-ymmsl>=0.12.0,<0.13
+git+https://github.com/multiscale/ymmsl-python.git@develop#egg=ymmsl
 qcg-pilotjob==0.13.1
 
diff --git a/libmuscle/cpp/src/libmuscle/communicator.cpp b/libmuscle/cpp/src/libmuscle/communicator.cpp
index 6b522259..644d67ae 100644
--- a/libmuscle/cpp/src/libmuscle/communicator.cpp
+++ b/libmuscle/cpp/src/libmuscle/communicator.cpp
@@ -116,7 +116,7 @@ void Communicator::send_message(
         // log sending on disconnected port
         return;
 
-    // Port const & port = ports_.at(port_name);
+    Port & port = ports_.at(port_name);
 
     // TODO start profile event
 
@@ -126,14 +126,15 @@ void Communicator::send_message(
     Data settings_overlay(message.settings());
 
     Optional<int> port_length;
-    if (ports_.at(port_name).is_resizable())
-        port_length = ports_.at(port_name).get_length();
+    if (port.is_resizable())
+        port_length = port.get_length();
 
     for (auto recv_endpoint : recv_endpoints) {
         MPPMessage mpp_message(
                 snd_endpoint.ref(), recv_endpoint.ref(),
                 port_length, message.timestamp(), Optional<double>(),
-                settings_overlay, message.data());
+                settings_overlay, port.get_num_messages(slot), -1.0,
+                message.data());
 
         if (message.has_next_timestamp())
             mpp_message.next_timestamp = message.next_timestamp();
@@ -141,6 +142,9 @@ void Communicator::send_message(
         auto message_bytes = std::make_unique<DataConstRef>(mpp_message.encoded());
         post_office_.deposit(recv_endpoint.ref(), std::move(message_bytes));
     }
+
+    port.increment_num_messages(slot);
+
     // TODO: stop and complete profile event
 }
 
@@ -207,6 +211,33 @@ Message Communicator::receive_message(
 
     // TODO stop and finalise profile event
 
+    int expected_message_number = port.get_num_messages(slot);
+    // TODO: handle f_init port counts for STATELESS and WEAKLY_STATEFUL
+    // components which didn't load a snapshot
+    if (expected_message_number != mpp_message.message_number) {
+        if (expected_message_number - 1 == mpp_message.message_number and
+                port.is_resuming(slot)) {
+            if (slot.is_set())
+                logger_.debug("Discarding received message on ", port_name,
+                              "[", slot.get(), "]: resuming from weakly",
+                              " consistent snapshot");
+            else
+                logger_.debug("Discarding received message on ", port_name,
+                              ": resuming from weakly constistent snapshot");
+            port.set_resumed(slot);
+            return receive_message(port_name, slot, default_msg);
+        }
+        std::ostringstream oss;
+        oss << "Received message on " << port_name;
+        if (slot.is_set())
+            oss << "[" << slot.get() << "]";
+        oss << " with unexpected message number " << mpp_message.message_number;
+        oss << ". Was expecting " << expected_message_number;
+        oss << ". Are you resuming from an inconsistent snapshot?";
+        throw std::runtime_error(oss.str());
+    }
+    port.increment_num_messages(slot);
+
     if (slot.is_set())
         logger_.debug("Received message on ", port_name, "[", slot.get(), "]");
     else
diff --git a/libmuscle/cpp/src/libmuscle/mmp_client.cpp b/libmuscle/cpp/src/libmuscle/mmp_client.cpp
index acd6672d..de50e894 100644
--- a/libmuscle/cpp/src/libmuscle/mmp_client.cpp
+++ b/libmuscle/cpp/src/libmuscle/mmp_client.cpp
@@ -105,7 +105,7 @@ void MMPClient::register_instance(
 
     auto response = call_manager_(request);
 
-    if (response.size() > 1)
+    if (response[0].as<int>() == static_cast<int>(ResponseType::error))
         throw std::runtime_error(
                 "Error registering instance: " + response[1].as<std::string>());
 }
diff --git a/libmuscle/cpp/src/libmuscle/mpp_message.cpp b/libmuscle/cpp/src/libmuscle/mpp_message.cpp
index 2962e31c..344adb23 100644
--- a/libmuscle/cpp/src/libmuscle/mpp_message.cpp
+++ b/libmuscle/cpp/src/libmuscle/mpp_message.cpp
@@ -16,6 +16,7 @@ MPPMessage::MPPMessage(
             ::libmuscle::impl::Optional<int> port_length,
             double timestamp, ::libmuscle::impl::Optional<double> next_timestamp,
             DataConstRef const & settings_overlay,
+            int message_number, double saved_until,
             DataConstRef const & data
             )
         : sender(sender)
@@ -24,6 +25,8 @@ MPPMessage::MPPMessage(
         , timestamp(timestamp)
         , next_timestamp(next_timestamp)
         , settings_overlay(settings_overlay)
+        , message_number(message_number)
+        , saved_until(saved_until)
         , data(data)
     {}
 
@@ -48,6 +51,8 @@ MPPMessage MPPMessage::from_bytes(DataConstRef const & data) {
             dict["timestamp"].as<double>(),
             next_timestamp,
             dict["settings_overlay"],
+            dict["message_number"].as<int>(),
+            dict["saved_until"].as<double>(),
             dict["data"]);
 }
 
@@ -67,6 +72,8 @@ DataConstRef MPPMessage::encoded() const {
             "timestamp", timestamp,
             "next_timestamp", next_timestamp_data,
             "settings_overlay", settings_overlay,
+            "message_number", message_number,
+            "saved_until", saved_until,
             "data", data
             );
 
diff --git a/libmuscle/cpp/src/libmuscle/mpp_message.hpp b/libmuscle/cpp/src/libmuscle/mpp_message.hpp
index 50e8a49b..69f15645 100644
--- a/libmuscle/cpp/src/libmuscle/mpp_message.hpp
+++ b/libmuscle/cpp/src/libmuscle/mpp_message.hpp
@@ -32,8 +32,8 @@ struct MPPMessage {
             ::ymmsl::Reference const & sender, ::ymmsl::Reference const & receiver,
             ::libmuscle::impl::Optional<int> port_length,
             double timestamp, ::libmuscle::impl::Optional<double> next_timestamp,
-            DataConstRef const & settings_overlay,
-            DataConstRef const & data);
+            DataConstRef const & settings_overlay, int message_number,
+            double saved_until, DataConstRef const & data);
 
     /** Create an MCP Message from an encoded buffer.
      *
@@ -53,6 +53,8 @@ struct MPPMessage {
     double timestamp;
     ::libmuscle::impl::Optional<double> next_timestamp;
     DataConstRef settings_overlay;
+    int message_number;
+    double saved_until;
     DataConstRef data;
 };
 
diff --git a/libmuscle/cpp/src/libmuscle/port.cpp b/libmuscle/cpp/src/libmuscle/port.cpp
index 691c75cf..1dd32a44 100644
--- a/libmuscle/cpp/src/libmuscle/port.cpp
+++ b/libmuscle/cpp/src/libmuscle/port.cpp
@@ -11,6 +11,19 @@ using ymmsl::Identifier;
 using ymmsl::Operator;
 
 
+namespace {
+
+template< typename T>
+inline void extend_vector_to_size(
+        std::vector<T> &vec, const int minsize, const T &val) {
+    if(static_cast<int>(vec.size()) < minsize) {
+        vec.resize(minsize, val);
+    }
+}
+
+}
+
+
 namespace libmuscle { namespace impl {
 
 Port::Port(
@@ -53,6 +66,8 @@ Port::Port(
     }
 
     is_resizable_ = is_vector && (our_ndims == static_cast<int>(peer_dims.size()));
+    num_messages_.resize(std::max(1, length_), 0);
+    is_resuming_.resize(std::max(1, length_), false);
 }
 
 bool Port::is_connected() const {
@@ -94,6 +109,11 @@ void Port::set_length(int length) {
     if (length != length_) {
         length_ = length;
         is_open_ = std::vector<bool>(length_, true);
+        // Using extend here to not discard any information about message
+        // numbers between resizes. Note that _num_messages and _is_resuming
+        // may be longer than self._length!
+        extend_vector_to_size(num_messages_, std::max(1, length_), 0);
+        extend_vector_to_size(is_resuming_, std::max(1, length_), false);
     }
 }
 
@@ -105,5 +125,38 @@ void Port::set_closed(int slot) {
     is_open_[slot] = false;
 }
 
+void Port::restore_message_counts(const std::vector<int> &num_messages) {
+    num_messages_ = std::vector<int>(num_messages);
+    is_resuming_.clear();
+    is_resuming_.resize(num_messages_.size(), true);
+    extend_vector_to_size(num_messages_, std::max(1, length_), 0);
+    extend_vector_to_size(is_resuming_, std::max(1, length_), false);
+}
+
+const std::vector<int> & Port::get_message_counts() const {
+    return num_messages_;
+}
+
+void Port::increment_num_messages(Optional<int> slot) {
+    int s = slot.is_set() ? slot.get() : 0;
+    num_messages_[s] ++;
+    set_resumed(s);
+}
+
+int Port::get_num_messages(Optional<int> slot) const {
+    int s = slot.is_set() ? slot.get() : 0;
+    return num_messages_[s];
+}
+
+bool Port::is_resuming(Optional<int> slot) const {
+    int s = slot.is_set() ? slot.get() : 0;
+    return is_resuming_[s];
+}
+
+void Port::set_resumed(Optional<int> slot) {
+    int s = slot.is_set() ? slot.get() : 0;
+    is_resuming_[s] = false;
+}
+
 } }
 
diff --git a/libmuscle/cpp/src/libmuscle/port.hpp b/libmuscle/cpp/src/libmuscle/port.hpp
index 41bf8898..908fb270 100644
--- a/libmuscle/cpp/src/libmuscle/port.hpp
+++ b/libmuscle/cpp/src/libmuscle/port.hpp
@@ -11,6 +11,10 @@ namespace libmuscle { namespace impl {
  * Ports can be used to send or receive messages. They have a name and an
  * operator, as well as a set of dimensions that determines the valid slot
  * indices for sending or receiving on this port.
+ *
+ * Ports keep track of the amount of messages sent or received on the port.
+ * However, the actual incrementing and validation is done in
+ * Communicator.
  */
 class Port : public ::ymmsl::Port {
     public:
@@ -104,11 +108,50 @@ class Port : public ::ymmsl::Port {
          */
         void set_closed(int slot);
 
+        /** Restore message counts from a snapshot.
+         *
+         * @param num_messages message counts of the snapshot
+         */
+        void restore_message_counts(const std::vector<int> &num_messages);
+
+        /** Get the message counts for all slots in this port
+         */
+        const std::vector<int> & get_message_counts() const;
+
+        /** Increment amount of messages sent or received.
+         *
+         * @param slot The slot that is sent/received on
+         */
+        void increment_num_messages(Optional<int> slot = {});
+
+        /** Get the amount of messages sent or received
+         *
+         * @param slot The slot that is sent/received on
+         */
+        int get_num_messages(Optional<int> slot = {}) const;
+
+        /** True when this port has resumed.
+         *
+         * After resumption, each port/slot may discard exactly one message.
+         * is_resuming keeps track of this state.
+         *
+         * @param slot The slot that is sent/received on
+         */
+        bool is_resuming(Optional<int> slot = {}) const;
+
+        /** Mark that this port has resumed and may no longer discard messages.
+         *
+         * @param slot The slot that is sent/received on
+         */
+        void set_resumed(Optional<int> slot = {});
+
     private:
         bool is_connected_;
         int length_;
         bool is_resizable_;
         std::vector<bool> is_open_;
+        std::vector<int> num_messages_;
+        std::vector<bool> is_resuming_;
 };
 
 } }
diff --git a/libmuscle/python/libmuscle/api_guard.py b/libmuscle/python/libmuscle/api_guard.py
new file mode 100644
index 00000000..56012448
--- /dev/null
+++ b/libmuscle/python/libmuscle/api_guard.py
@@ -0,0 +1,248 @@
+from enum import auto, Enum
+from typing import Optional
+
+
+class APIPhase(Enum):
+    """Different phases that the user code traverses.
+
+    These values describe different regions that the model code can be
+    in for the case where checkpointing is implemented. By tracking
+    the phase that the model should be in, we can detect incorrect API
+    usage.
+
+    This does not match the yMMSL operators, as it is more
+    fine-grained and concerns checkpointing, which is not represented
+    in the SEL.
+
+    Note that AFTER_REUSE_INSTANCE and BEFORE_RESUMING refer to the
+    same place in the code. AFTER_REUSE_INSTANCE is used when we
+    don't know yet if the code has checkpointing support, and so we
+    don't know whether the next call is to resuming() or to
+    reuse_instance(). Once a checkpointing function has been called,
+    we know that we should expect resume() after reuse_instance() and
+    we use BEFORE_RESUMING accordingly.
+    """
+    BEFORE_FIRST_REUSE_INSTANCE = auto()
+    """Before the first time calling reuse_instance"""
+
+    BEFORE_REUSE_INSTANCE = auto()
+    """Before calling reuse_instance"""
+
+    AFTER_REUSE_INSTANCE = auto()
+    """At the top of the reuse loop"""
+
+    BEFORE_RESUMING = auto()
+    """Between reuse_instance and resuming"""
+
+    BEFORE_LOAD_SNAPSHOT = auto()
+    """Between resuming and load_snapshot"""
+
+    BEFORE_SHOULD_INIT = auto()
+    """After resuming, before should_init"""
+
+    BEFORE_SHOULD_SAVE_SNAPSHOT = auto()
+    """Between should_init and should_save*"""
+
+    BEFORE_SAVE_SNAPSHOT = auto()
+    """Between should_save_snapshot and save_snapshot"""
+
+    BEFORE_SAVE_FINAL_SNAPSHOT = auto()
+    """Between should_save_final_snapshot and save_final_snapshot"""
+
+    AFTER_REUSE_LOOP = auto()
+    """After the final call to reuse_instance()"""
+
+
+class APIGuard:
+    """Keeps track of and checks in which phase the model is.
+
+    The verify_* functions are called when the corresponding function
+    on Instance is called, to check that we're in the right phase. They
+    raise a RuntimeError if there's a problem. The *_done functions are
+    called to signal that the corresponding function finished
+    successfully, and that we are moving on to the next phase.
+    """
+    def __init__(self) -> None:
+        """Create an APIPhaseTracker.
+
+        This starts the tracker in BEFORE_FIRST_REUSE_INSTANCE.
+        """
+        self._phase = APIPhase.BEFORE_FIRST_REUSE_INSTANCE
+        self._uses_checkpointing = None     # type: Optional[bool]
+
+    def _generic_error_messages(self, verify_phase: str) -> None:
+        if self._phase in (
+                APIPhase.BEFORE_FIRST_REUSE_INSTANCE,
+                APIPhase.AFTER_REUSE_LOOP):
+            msg = f'Please only call {verify_phase} inside the reuse loop.'
+        elif self._phase == APIPhase.BEFORE_REUSE_INSTANCE:
+            msg = (
+                    'Please do not call {verify_phase} after'
+                    ' should_save_final_snapshot. should_save_final_snapshot'
+                    ' should be at the end of the reuse loop.')
+        elif self._phase == APIPhase.AFTER_REUSE_INSTANCE:
+            msg = (
+                    'Please call resuming first in the reuse loop, before'
+                    f' {verify_phase}')
+        elif self._phase == APIPhase.BEFORE_RESUMING:
+            msg = 'Inside the reuse loop you must call resuming first.'
+        elif self._phase == APIPhase.BEFORE_LOAD_SNAPSHOT:
+            msg = (
+                    'If resuming returns True, then you must call'
+                    ' load_snapshot first.')
+        elif self._phase == APIPhase.BEFORE_SHOULD_INIT:
+            msg = 'After calling resuming, you must call should_init first.'
+        elif self._phase == APIPhase.BEFORE_SHOULD_SAVE_SNAPSHOT:
+            msg = 'You must call save_snapshot or save_final_snapshot first.'
+        elif self._phase == APIPhase.BEFORE_SAVE_SNAPSHOT:
+            msg = (
+                    'If should_save_snapshot returns True, then you must'
+                    ' call save_snapshot first.')
+        elif self._phase == APIPhase.BEFORE_SAVE_FINAL_SNAPSHOT:
+            msg = (
+                    'If should_save_final_snapshot returns True, then you'
+                    ' must call save_final_snapshot first.')
+        else:
+            return
+        raise RuntimeError(msg)
+
+    def uses_checkpointing(self) -> bool:
+        """Return whether the code is using checkpointing.
+
+        We can only determine that the code doesn't use checkpointing
+        if there are no checkpointing calls between the first and
+        second calls to reuse_instance. So this function should only
+        be called after the second call to verify_reuse_instance, or
+        it may raise if the code does not use checkpointing.
+
+        Raises:
+            RuntimeError: if we are at a point where we cannot know
+                the answer yet.
+        """
+        if self._uses_checkpointing is not None:
+            return self._uses_checkpointing
+        raise RuntimeError(
+                'The API was implemented incorrectly, please consult the'
+                ' documentation.')
+
+    def verify_reuse_instance(self) -> None:
+        """Check reuse_instance()"""
+        if self._phase == APIPhase.AFTER_REUSE_INSTANCE:
+            self._uses_checkpointing = False
+        elif self._phase not in (
+                APIPhase.BEFORE_REUSE_INSTANCE,
+                APIPhase.BEFORE_FIRST_REUSE_INSTANCE):
+            raise RuntimeError(
+                    'We reached the end of the reuse loop without checking'
+                    ' if a snapshot should be saved. Please add at least'
+                    ' a should_save_final_snapshot and save_final_snapshot.')
+
+    def reuse_instance_done(self, reusing: bool) -> None:
+        """Update phase on successful reuse_instance().
+
+        Args:
+            reusing: Whether we are reusing or not.
+        """
+        if not reusing:
+            self._phase = APIPhase.AFTER_REUSE_LOOP
+        else:
+            if self._uses_checkpointing is None:
+                self._phase = APIPhase.AFTER_REUSE_INSTANCE
+            elif self._uses_checkpointing:
+                self._phase = APIPhase.BEFORE_RESUMING
+            else:
+                self._phase = APIPhase.BEFORE_REUSE_INSTANCE
+
+    def verify_resuming(self) -> None:
+        """Check resuming()"""
+        if self._phase not in (
+                APIPhase.BEFORE_RESUMING, APIPhase.AFTER_REUSE_INSTANCE):
+            raise RuntimeError(
+                    'Please call resuming() only as the first thing in the'
+                    ' reuse loop.')
+
+    def resuming_done(self, resuming: bool) -> None:
+        """Update phase on successful resuming().
+
+        Args:
+            resuming: Whether we're resuming or not.
+        """
+        self._uses_checkpointing = True
+        if resuming:
+            self._phase = APIPhase.BEFORE_LOAD_SNAPSHOT
+        else:
+            self._phase = APIPhase.BEFORE_SHOULD_INIT
+
+    def verify_load_snapshot(self) -> None:
+        """Check load_snapshot()"""
+        if self._phase != APIPhase.BEFORE_LOAD_SNAPSHOT:
+            raise RuntimeError(
+                    'Please check that we are resuming by calling resuming()'
+                    ' before calling load_snapshot()')
+
+    def load_snapshot_done(self) -> None:
+        """Update phase on successful load_snapshot()"""
+        self._phase = APIPhase.BEFORE_SHOULD_INIT
+
+    def verify_should_init(self) -> None:
+        """Check should_init()"""
+        if self._phase != APIPhase.BEFORE_SHOULD_INIT:
+            raise RuntimeError(
+                    'Please check whether to run f_init using should_init()'
+                    ' after resuming, and before trying to save a snapshot.')
+
+    def should_init_done(self) -> None:
+        """Update phase on successful should_init()"""
+        self._phase = APIPhase.BEFORE_SHOULD_SAVE_SNAPSHOT
+
+    def verify_should_save_snapshot(self) -> None:
+        """Check should_save_snapshot()"""
+        if self._phase != APIPhase.BEFORE_SHOULD_SAVE_SNAPSHOT:
+            self._generic_error_messages('should_save_snapshot')
+            raise RuntimeError()  # should be unreachable
+
+    def should_save_snapshot_done(self, should_save: bool) -> None:
+        """Update phase on successful should_save_snapshot().
+
+        Args:
+            should_save: Whether we should save or not.
+        """
+        if should_save:
+            self._phase = APIPhase.BEFORE_SAVE_SNAPSHOT
+
+    def verify_save_snapshot(self) -> None:
+        """Check save_snapshot()"""
+        if self._phase != APIPhase.BEFORE_SAVE_SNAPSHOT:
+            self._generic_error_messages('save_snapshot')
+            raise RuntimeError()  # should be unreachable
+
+    def save_snapshot_done(self) -> None:
+        """Update phase on successful save_snapshot()"""
+        self._phase = APIPhase.BEFORE_SHOULD_SAVE_SNAPSHOT
+
+    def verify_should_save_final_snapshot(self) -> None:
+        """Check should_save_final_snapshot()."""
+        if self._phase != APIPhase.BEFORE_SHOULD_SAVE_SNAPSHOT:
+            self._generic_error_messages('should_save_final_snapshot')
+            raise RuntimeError()  # should be unreachable
+
+    def should_save_final_snapshot_done(self, should_save: bool) -> None:
+        """Update phase on successful should_save_snapshot().
+
+        Args:
+            should_save: Whether we should save or not.
+        """
+        if should_save:
+            self._phase = APIPhase.BEFORE_SAVE_FINAL_SNAPSHOT
+        else:
+            self._phase = APIPhase.BEFORE_REUSE_INSTANCE
+
+    def verify_save_final_snapshot(self) -> None:
+        """Check should_save_final_snapshot()"""
+        if self._phase != APIPhase.BEFORE_SAVE_FINAL_SNAPSHOT:
+            self._generic_error_messages('save_final_snapshot')
+            raise RuntimeError()  # should be unreachable
+
+    def save_final_snapshot_done(self) -> None:
+        """Updates state on successful save_final_snapshot()"""
+        self._phase = APIPhase.BEFORE_REUSE_INSTANCE
diff --git a/libmuscle/python/libmuscle/checkpoint_triggers.py b/libmuscle/python/libmuscle/checkpoint_triggers.py
new file mode 100644
index 00000000..b0064c75
--- /dev/null
+++ b/libmuscle/python/libmuscle/checkpoint_triggers.py
@@ -0,0 +1,302 @@
+import bisect
+import logging
+import time
+from typing import List, Optional, Union
+
+from ymmsl import (
+        CheckpointRangeRule, CheckpointAtRule, CheckpointRule, Checkpoints)
+
+
+_logger = logging.getLogger(__name__)
+
+
+class CheckpointTrigger:
+    """Represents a trigger for creating snapshots"""
+
+    def next_checkpoint(self, cur_time: float) -> Optional[float]:
+        """Calculate the next checkpoint time
+
+        Args:
+            cur_time: current time.
+
+        Returns:
+            The time when a next checkpoint should be taken, or None if this
+            trigger has no checkpoint after cur_time.
+        """
+        raise NotImplementedError()
+
+    def previous_checkpoint(self, cur_time: float) -> Optional[float]:
+        """Calculate the previous checkpoint time
+
+        Args:
+            cur_time: current time.
+
+        Returns:
+            The time when a previous checkpoint should have been taken, or None
+            if this trigger has no checkpoint after cur_time.
+        """
+        raise NotImplementedError()
+
+
+class AtCheckpointTrigger(CheckpointTrigger):
+    """Represents a trigger based on an "at" checkpoint rule
+
+    This triggers at the specified times.
+    """
+
+    def __init__(self, at_rules: List[CheckpointAtRule]) -> None:
+        """Create an "at" checkpoint trigger
+
+        Args:
+            at: list of checkpoint moments
+        """
+        self._at = sorted([a for r in at_rules for a in r.at])
+
+    def next_checkpoint(self, cur_time: float) -> Optional[float]:
+        if cur_time >= self._at[-1]:
+            return None  # no future checkpoint left
+        idx = bisect.bisect(self._at, cur_time)
+        return self._at[idx]
+
+    def previous_checkpoint(self, cur_time: float) -> Optional[float]:
+        if cur_time < self._at[0]:
+            return None  # no previous checkpoint
+        idx = bisect.bisect(self._at, cur_time)
+        return self._at[idx - 1]
+
+
+class RangeCheckpointTrigger(CheckpointTrigger):
+    """Represents a trigger based on a "ranges" checkpoint rule
+
+    This triggers at a range of checkpoint moments.
+
+    Equivalent an "at" rule ``[start, start + step, start + 2*step, ...]`` for
+    as long as ``start + i*step <= stop``.
+
+    Stop may be omitted, in which case the range is infinite.
+
+    Start may be omitted, in which case the range is equivalent to an "at" rule
+    ``[..., -n*step, ..., -step, 0, step, 2*step, ...]`` for as long as
+    ``i*step <= stop``.
+
+    Note: the "every" rule is a special case of a range with start and stop
+    omitted, and is handled by this class as well
+    """
+
+    def __init__(self, range: CheckpointRangeRule) -> None:
+        """Create a range of checkpoints
+
+        Args:
+            range: checkpoint range defining start, stop and step.
+        """
+        self._start = range.start
+        self._stop = range.stop
+        self._every = range.every
+        self._last = None  # type: Union[int, float, None]
+        if self._stop is not None:
+            start = 0 if self._start is None else self._start
+            diff = self._stop - start
+            self._last = start + (diff // self._every) * self._every
+
+    def next_checkpoint(self, cur_time: float) -> Optional[float]:
+        if self._start is not None and cur_time < self._start:
+            return float(self._start)
+        if self._last is not None and cur_time >= self._last:
+            return None
+        start = 0 if self._start is None else self._start
+        diff = cur_time - start
+        return float(start + (diff // self._every + 1) * self._every)
+
+    def previous_checkpoint(self, cur_time: float) -> Optional[float]:
+        if self._start is not None and cur_time < self._start:
+            return None
+        if self._last is not None and cur_time > self._last:
+            return float(self._last)
+        start = 0 if self._start is None else self._start
+        diff = cur_time - start
+        return float(start + (diff // self._every) * self._every)
+
+
+class CombinedCheckpointTriggers(CheckpointTrigger):
+    """Checkpoint trigger based on a combination of "every", "at" and "ranges"
+    """
+
+    def __init__(self, checkpoint_rules: List[CheckpointRule]) -> None:
+        """Create a new combined checkpoint trigger from the given rules
+
+        Args:
+            checkpoint_rules: checkpoint rules (from ymmsl)
+        """
+        self._triggers = []     # type: List[CheckpointTrigger]
+        at_rules = []           # type: List[CheckpointAtRule]
+        for rule in checkpoint_rules:
+            if isinstance(rule, CheckpointAtRule):
+                if rule.at:
+                    at_rules.append(rule)
+            elif isinstance(rule, CheckpointRangeRule):
+                self._triggers.append(RangeCheckpointTrigger(rule))
+            else:
+                raise RuntimeError('Unknown checkpoint rule')
+        if at_rules:
+            self._triggers.append(AtCheckpointTrigger(at_rules))
+
+    def next_checkpoint(self, cur_time: float) -> Optional[float]:
+        checkpoints = (trigger.next_checkpoint(cur_time)
+                       for trigger in self._triggers)
+        # return earliest of all not-None next-checkpoints
+        return min((checkpoint
+                    for checkpoint in checkpoints
+                    if checkpoint is not None),
+                   default=None)  # return None if all triggers return None
+
+    def previous_checkpoint(self, cur_time: float) -> Optional[float]:
+        checkpoints = (trigger.previous_checkpoint(cur_time)
+                       for trigger in self._triggers)
+        # return latest of all not-None previous-checkpoints
+        return max((checkpoint
+                    for checkpoint in checkpoints
+                    if checkpoint is not None),
+                   default=None)  # return None if all triggers return None
+
+
+class TriggerManager:
+    """Manages all checkpoint triggers and checks if a snapshot must be saved.
+    """
+
+    def __init__(self) -> None:
+        self._has_checkpoints = False
+        self._last_triggers = []    # type: List[str]
+        self._cpts_considered_until = float('-inf')
+
+    def set_checkpoint_info(
+            self, elapsed: float, checkpoints: Checkpoints) -> None:
+        """Register checkpoint info received from the muscle manager.
+        """
+        self._mono_to_elapsed = elapsed - time.monotonic()
+
+        if not checkpoints:
+            self._has_checkpoints = False
+            return
+
+        self._has_checkpoints = True
+
+        self._checkpoint_at_end = checkpoints.at_end
+
+        self._wall = CombinedCheckpointTriggers(checkpoints.wallclock_time)
+        self._prevwall = 0.0
+        self._nextwall = self._wall.next_checkpoint(0.0)  # type: Optional[float]
+
+        self._sim = CombinedCheckpointTriggers(checkpoints.simulation_time)
+        self._prevsim = None        # type: Optional[float]
+        self._nextsim = None        # type: Optional[float]
+        self._sim_reset = True
+
+    def elapsed_walltime(self) -> float:
+        """Returns elapsed wallclock_time in seconds.
+        """
+        return time.monotonic() + self._mono_to_elapsed
+
+    def checkpoints_considered_until(self) -> float:
+        """Return elapsed time of last should_save*
+        """
+        return self._cpts_considered_until
+
+    def harmonise_wall_time(self, at_least: float) -> None:
+        """Ensure our elapsed time is at least the given value
+        """
+        cur = self.elapsed_walltime()
+        if cur < at_least:
+            _logger.debug(
+                    'Harmonise wall time: advancing clock by %f seconds',
+                    at_least - cur)
+            self._mono_to_elapsed += at_least - cur
+
+    def snapshots_enabled(self) -> bool:
+        """Check if the current workflow has snapshots enabled.
+        """
+        return self._has_checkpoints
+
+    def should_save_snapshot(self, timestamp: float) -> bool:
+        """Handles instance.should_save_snapshot
+        """
+        if not self._has_checkpoints:
+            return False
+
+        return self.__should_save(timestamp)
+
+    def should_save_final_snapshot(
+            self, do_reuse: bool, f_init_max_timestamp: Optional[float]
+            ) -> bool:
+        """Handles instance.should_save_final_snapshot
+        """
+        if not self._has_checkpoints:
+            return False
+
+        value = False
+        if not do_reuse:
+            if self._checkpoint_at_end:
+                value = True
+                self._last_triggers.append('at_end')
+        elif f_init_max_timestamp is None:
+            # No F_INIT messages received: reuse triggered on muscle_settings_in
+            # message.
+            _logger.debug('Reuse triggered by muscle_settings_in.'
+                          ' Not creating a snapshot.')
+            self._sim_reset = True
+        else:
+            value = self.__should_save(f_init_max_timestamp)
+
+        return value
+
+    def update_checkpoints(self, timestamp: float) -> None:
+        """Update last and next checkpoint times when a snapshot is made.
+
+        Args:
+            timestamp: timestamp as reported by the instance (or from incoming
+                F_INIT messages for save_final_snapshot).
+        """
+        self._prevwall = self.elapsed_walltime()
+        self._nextwall = self._wall.next_checkpoint(self._prevwall)
+
+        self._prevsim = timestamp
+        self._nextsim = self._sim.next_checkpoint(timestamp)
+
+        # this method is also called during resume, after which we no longer
+        # consider the simulation_time as reset
+        self._sim_reset = False
+
+    def get_triggers(self) -> List[str]:
+        """Get trigger description(s) for the current reason for checkpointing.
+        """
+        triggers = self._last_triggers
+        self._last_triggers = []
+        return triggers
+
+    def __should_save(self, simulation_time: float) -> bool:
+        """Check if a checkpoint should be taken
+
+        Args:
+            simulation_time: current/next timestamp as reported by the instance
+        """
+        if self._sim_reset:
+            # we cannot make assumptions about the start time of a simulation,
+            # a t=-1000 could make sense if t represents years since CE
+            # and we should not disallow checkpointing for negative t
+            previous = self._sim.previous_checkpoint(simulation_time)
+            if previous is not None:
+                # there is a checkpoint rule before the current moment, assume
+                # we should have taken a snapshot back then
+                self._nextsim = previous
+            else:
+                self._nextsim = self._sim.next_checkpoint(simulation_time)
+            self._sim_reset = False
+
+        walltime = self.elapsed_walltime()
+        self._cpts_considered_until = walltime
+
+        self._last_triggers = []
+        if self._nextwall is not None and walltime >= self._nextwall:
+            self._last_triggers.append(f"wallclock_time >= {self._nextwall}")
+        if self._nextsim is not None and simulation_time >= self._nextsim:
+            self._last_triggers.append(f"simulation_time >= {self._nextsim}")
+        return bool(self._last_triggers)
diff --git a/libmuscle/python/libmuscle/communicator.py b/libmuscle/python/libmuscle/communicator.py
index 46b76468..69272f78 100644
--- a/libmuscle/python/libmuscle/communicator.py
+++ b/libmuscle/python/libmuscle/communicator.py
@@ -181,7 +181,8 @@ def get_port(self, port_name: str) -> Port:
 
     def send_message(
             self, port_name: str, message: Message,
-            slot: Optional[int] = None) -> None:
+            slot: Optional[int] = None,
+            checkpoints_considered_until: float = float('-inf')) -> None:
         """Send a message and settings to the outside world.
 
         Sending is non-blocking, a copy of the message will be made
@@ -191,6 +192,8 @@ def send_message(
             port_name: The port on which this message is to be sent.
             message: The message to be sent.
             slot: The slot to send the message on, if any.
+            checkpoints_considered_until: When we last checked if we
+                should save a snapshot (wallclock time).
         """
         if slot is None:
             _logger.debug('Sending message on {}'.format(port_name))
@@ -218,18 +221,22 @@ def send_message(
                 snd_endpoint.port, slot_list)
 
         port_length = None
-        if self._ports[port_name].is_resizable():
-            port_length = self._ports[port_name].get_length()
+        if port.is_resizable():
+            port_length = port.get_length()
 
         for recv_endpoint in recv_endpoints:
             mcp_message = MPPMessage(snd_endpoint.ref(), recv_endpoint.ref(),
                                      port_length,
                                      message.timestamp, message.next_timestamp,
                                      cast(Settings, message.settings),
+                                     port.get_num_messages(slot),
+                                     checkpoints_considered_until,
                                      message.data)
             encoded_message = mcp_message.encoded()
             self._post_office.deposit(recv_endpoint.ref(), encoded_message)
 
+        port.increment_num_messages(slot)
+
         profile_event.stop()
         if port.is_vector():
             profile_event.port_length = port.get_length()
@@ -237,7 +244,7 @@ def send_message(
 
     def receive_message(self, port_name: str, slot: Optional[int] = None,
                         default: Optional[Message] = None
-                        ) -> Message:
+                        ) -> Tuple[Message, float]:
         """Receive a message and attached settings overlay.
 
         Receiving is a blocking operation. This function will contact
@@ -257,19 +264,20 @@ def receive_message(self, port_name: str, slot: Optional[int] = None,
         Returns:
             The received message, with message.settings holding
             the settings overlay. The settings attribute is
-            guaranteed to not be None.
+            guaranteed to not be None. Secondly, the saved_until
+            metadata field from the received message.
 
         Raises:
             RuntimeError: If no default was given and the port is not
                 connected.
         """
         if slot is None:
-            _logger.debug('Waiting for message on {}'.format(port_name))
+            port_and_slot = port_name
             slot_list = []      # type: List[int]
         else:
-            _logger.debug('Waiting for message on {}[{}]'.format(
-                port_name, slot))
+            port_and_slot = f"{port_name}[{slot}]"
             slot_list = [slot]
+        _logger.debug('Waiting for message on {}'.format(port_and_slot))
 
         recv_endpoint = self.__get_endpoint(port_name, slot_list)
 
@@ -283,7 +291,7 @@ def receive_message(self, port_name: str, slot: Optional[int] = None,
             _logger.debug(
                     'No message received on {} as it is not connected'.format(
                         port_name))
-            return default
+            return default, float('-inf')
 
         if port_name in self._ports:
             port = self._ports[port_name]
@@ -301,36 +309,45 @@ def receive_message(self, port_name: str, slot: Optional[int] = None,
         snd_endpoint = self._peer_manager.get_peer_endpoints(
                 recv_endpoint.port, slot_list)[0]
         client = self.__get_client(snd_endpoint.instance())
-        mcp_message_bytes = client.receive(recv_endpoint.ref())
-        mcp_message = MPPMessage.from_bytes(mcp_message_bytes)
+        mpp_message_bytes = client.receive(recv_endpoint.ref())
+        mpp_message = MPPMessage.from_bytes(mpp_message_bytes)
 
-        if mcp_message.port_length is not None:
+        if mpp_message.port_length is not None:
             if port.is_resizable():
-                port.set_length(mcp_message.port_length)
+                port.set_length(mpp_message.port_length)
 
-        if isinstance(mcp_message.data, ClosePort):
+        if isinstance(mpp_message.data, ClosePort):
             port.set_closed(slot)
 
         message = Message(
-                mcp_message.timestamp, mcp_message.next_timestamp,
-                mcp_message.data, mcp_message.settings_overlay)
+                mpp_message.timestamp, mpp_message.next_timestamp,
+                mpp_message.data, mpp_message.settings_overlay)
 
         profile_event.stop()
         if port.is_vector():
             profile_event.port_length = port.get_length()
-        profile_event.message_size = len(mcp_message_bytes)
-
-        if slot is None:
-            _logger.debug('Received message on {}'.format(port_name))
-            if isinstance(mcp_message.data, ClosePort):
-                _logger.debug('Port {} is now closed'.format(port_name))
-        else:
-            _logger.debug('Received message on {}[{}]'.format(port_name, slot))
-            if isinstance(mcp_message.data, ClosePort):
-                _logger.debug('Port {}[{}] is now closed'.format(
-                    port_name, slot))
-
-        return message
+        profile_event.message_size = len(mpp_message_bytes)
+
+        expected_message_number = port.get_num_messages(slot)
+        if expected_message_number != mpp_message.message_number:
+            if (expected_message_number - 1 == mpp_message.message_number and
+                    port.is_resuming(slot)):
+                _logger.debug(f'Discarding received message on {port_and_slot}'
+                              ': resuming from weakly consistent snapshot')
+                port.set_resumed(slot)
+                return self.receive_message(port_name, slot, default)
+            raise RuntimeError(f'Received message on {port_and_slot} with'
+                               ' unexpected message number'
+                               f' {mpp_message.message_number}. Was expecting'
+                               f' {expected_message_number}. Are you resuming'
+                               ' from an inconsistent snapshot?')
+        port.increment_num_messages(slot)
+
+        _logger.debug('Received message on {}'.format(port_and_slot))
+        if isinstance(mpp_message.data, ClosePort):
+            _logger.debug('Port {} is now closed'.format(port_and_slot))
+
+        return message, mpp_message.saved_until
 
     def close_port(self, port_name: str, slot: Optional[int] = None
                    ) -> None:
@@ -361,6 +378,29 @@ def shutdown(self) -> None:
         for server in self._servers:
             server.close()
 
+    def restore_message_counts(self, port_message_counts: Dict[str, List[int]]
+                               ) -> None:
+        """Restore message counts on all ports
+        """
+        for port_name, num_messages in port_message_counts.items():
+            if port_name == "muscle_settings_in":
+                self._muscle_settings_in.restore_message_counts(num_messages)
+            elif port_name in self._ports:
+                self._ports[port_name].restore_message_counts(num_messages)
+            else:
+                raise RuntimeError(f'Unknown port {port_name} in snapshot.'
+                                   ' Have your port definitions changed since'
+                                   ' the snapshot was taken?')
+
+    def get_message_counts(self) -> Dict[str, List[int]]:
+        """Get message counts for all ports on the communicator
+        """
+        port_message_counts = {port_name: port.get_message_counts()
+                               for port_name, port in self._ports.items()}
+        port_message_counts["muscle_settings_in"] = \
+            self._muscle_settings_in.get_message_counts()
+        return port_message_counts
+
     def __instance_id(self) -> Reference:
         """Returns our complete instance id.
         """
diff --git a/libmuscle/python/libmuscle/grid.py b/libmuscle/python/libmuscle/grid.py
index 34efc806..bfe26715 100644
--- a/libmuscle/python/libmuscle/grid.py
+++ b/libmuscle/python/libmuscle/grid.py
@@ -6,13 +6,13 @@
 class Grid:
     """Represents a grid of data to send or receive.
 
-    Note that for received grids, the array of data is a read-only
-    NumPy array. If you have another array that you want to put the
-    received data into, use ``np.copyto(dest, source)`` to copy the
-    contents of the received array across into your destination array.
-    If you don't have an array yet and want a writable version of the
-    received array, use ``array.copy()`` to create a writable copy.
-    See the tutorial for examples.
+    Note that for received grids, the array of data is a read-only NumPy array.
+    If you have another array that you want to put the received data into, use
+    :external:py:func:`np.copyto(dest, source) <numpy.copyto>` to copy the
+    contents of the received array across into your destination array. If you
+    don't have an array yet and want a writable version of the received array,
+    use :external:py:meth:`array.copy()<numpy.ndarray.copy>` to create a
+    writable copy. See the tutorial for examples.
 
     Attributes:
         array (np.ndarray): An array of data
@@ -26,10 +26,12 @@ def __init__(
         A Grid object represents an multi-dimensional array of data. It
         has a type, a shape, and optionally a list of index names.
 
-        Supported data types are 4- and 8-byte integers (numpy.int32,
-        numpy.int64), 4- and 8-byte floats (numpy.float32,
-        numpy.float64), and booleans (np.bool_, np.bool8). The ``data``
-        argument must be a NumPy array of one of those types.
+        Supported data types are 4- and 8-byte integers
+        (:external:py:attr:`numpy.int32`, :external:py:attr:`numpy.int64`),
+        4- and 8-byte floats (:external:py:attr:`numpy.float32`,
+        :external:py:attr:`numpy.float64`), and booleans
+        (:external:py:class:`numpy.bool_`, :external:py:attr:`numpy.bool8`). The
+        ``data`` argument must be a NumPy array of one of those types.
 
         If ``indexes`` is given, then it must be a list of strings of
         the same length as the number of dimensions of ``data``, and
diff --git a/libmuscle/python/libmuscle/instance.py b/libmuscle/python/libmuscle/instance.py
index ae952111..e4e685f7 100644
--- a/libmuscle/python/libmuscle/instance.py
+++ b/libmuscle/python/libmuscle/instance.py
@@ -7,8 +7,10 @@
 from typing_extensions import Literal
 
 from ymmsl import (Identifier, Operator, SettingValue, Port, Reference,
-                   Settings)
+                   Settings, KeepsStateForNextUse)
 
+from libmuscle.api_guard import APIGuard
+from libmuscle.checkpoint_triggers import TriggerManager
 from libmuscle.communicator import Communicator, Message
 from libmuscle.settings_manager import SettingsManager
 from libmuscle.logging import LogLevel
@@ -17,6 +19,7 @@
 from libmuscle.mmp_client import MMPClient
 from libmuscle.profiler import Profiler
 from libmuscle.profiling import ProfileEventType
+from libmuscle.snapshot_manager import SnapshotManager
 from libmuscle.util import extract_log_file_location
 
 
@@ -32,16 +35,24 @@ class Instance:
     This class provides a low-level send/receive API for the instance
     to use.
     """
-    def __init__(self, ports: Optional[Dict[Operator, List[str]]] = None
-                 ) -> None:
+    def __init__(
+            self, ports: Optional[Dict[Operator, List[str]]] = None,
+            keeps_state_for_next_use: KeepsStateForNextUse
+            = KeepsStateForNextUse.NECESSARY) -> None:
         """Create an Instance.
 
         Args:
-            ports: A list of port names for each operator of this
-                component.
+            ports: A list of port names for each
+                :external:py:class:`~ymmsl.Operator` of this component.
+            keeps_state_for_next_use: Indicate whether this instance carries
+                state between iterations of the reuse loop. See
+                :external:py:class:`ymmsl.KeepsStateForNextUse` for a
+                description of the options.
         """
         self.__is_shut_down = False
 
+        self._keeps_state = KeepsStateForNextUse(keeps_state_for_next_use)
+
         # Note that these are accessed by Muscle3, but otherwise private.
         self._name, self._index = self.__make_full_name()
         """Name and index of this instance."""
@@ -52,6 +63,9 @@ def __init__(self, ports: Optional[Dict[Operator, List[str]]] = None
 
         self.__set_up_logging()
 
+        self._api_guard = APIGuard()
+        """Checks that the user uses the API correctly."""
+
         self._profiler = Profiler(self._instance_name(), self.__manager)
         """Profiler for this instance."""
 
@@ -65,13 +79,50 @@ def __init__(self, ports: Optional[Dict[Operator, List[str]]] = None
         self._settings_manager = SettingsManager()
         """Settings for this instance."""
 
-        self._first_run = True
-        """Keeps track of whether this is the first reuse run."""
+        self._snapshot_manager = SnapshotManager(
+                self._instance_name(), self.__manager, self._communicator)
+        """Resumes, loads and saves snapshots."""
+
+        self._trigger_manager = TriggerManager()
+        """Keeps track of checkpoints and triggers snapshots."""
+
+        self._first_run = None          # type: Optional[bool]
+        """Whether this is the first iteration of the reuse loop"""
+
+        self._do_reuse = None           # type: Optional[bool]
+        """Whether to enter this iteration of the reuse loop
+
+        This is None during the reuse loop, and set between
+        should_save_final_snapshot and reuse_instance.
+        """
+
+        self._do_resume = False
+        """Whether to resume on this iteration of the reuse loop"""
+
+        self._do_init = False
+        """Whether to do f_init on this iteration of the reuse loop"""
 
         self._f_init_cache = dict()     # type: _FInitCacheType
+        """Stores pre-received messages for f_init ports"""
 
         self._register()
         self._connect()
+
+        # Note: get_checkpoint_info needs to have the ports initialized
+        # so it comes after self._connect()
+        checkpoint_info = self.__manager.get_checkpoint_info(
+                self._instance_name())
+
+        elapsed_time, checkpoints = checkpoint_info[0:2]
+        self._trigger_manager.set_checkpoint_info(elapsed_time, checkpoints)
+
+        resume_snapshot, snapshot_dir = checkpoint_info[2:4]
+        saved_at = self._snapshot_manager.prepare_resume(
+                resume_snapshot, snapshot_dir)
+
+        if saved_at is not None:
+            self._trigger_manager.update_checkpoints(saved_at)
+
         self._set_local_log_level()
         self._set_remote_log_level()
 
@@ -103,40 +154,51 @@ def reuse_instance(self, apply_overlay: bool = True) -> bool:
                 overlay or to save it. If you're going to use
                 :meth:`receive_with_settings` on your F_INIT ports,
                 set this to False. If you don't know what that means,
-                just call `reuse_instance()` without specifying this
+                just call :meth:`reuse_instance()` without specifying this
                 and everything will be fine. If it turns out that you
                 did need to specify False, MUSCLE3 will tell you about
                 it in an error message and you can add it still.
+
+        Raises:
+            RuntimeError:
+                When implementing the checkpointing API, but libmuscle detected
+                incorrect API calls. The description of the RuntimeError
+                indicates which calls are incorrect or missing. For more
+                information see the checkpointing API documentation in
+                :meth:`resuming`, :meth:`load_snapshot`,
+                :meth:`should_save_snapshot`, :meth:`save_snapshot`,
+                :meth:`should_save_final_snapshot` and
+                :meth:`save_final_snapshot`, or the checkpointing tutorial.
         """
-        do_reuse = self.__receive_settings()
+        self._api_guard.verify_reuse_instance()
 
-        # TODO: _f_init_cache should be empty here, or the user didn't
-        # receive something that was sent on the last go-around.
-        # At least emit a warning.
-        self.__pre_receive_f_init(apply_overlay)
+        if self._do_reuse is not None:
+            # thank you, should_save_final_snapshot, for running this already
+            do_reuse = self._do_reuse
+            self._do_reuse = None
+        else:
+            do_reuse = self._decide_reuse_instance(apply_overlay)
 
-        self._set_local_log_level()
-        self._set_remote_log_level()
+        # now _first_run, _do_resume and _do_init are also set correctly
 
-        ports = self._communicator.list_ports()
-        f_init_not_connected = all(
-                [not self.is_connected(port)
-                 for port in ports.get(Operator.F_INIT, [])])
-        no_settings_in = not self._communicator.settings_in_connected()
+        do_implicit_checkpoint = (
+                not self._first_run and
+                not self._api_guard.uses_checkpointing() and
+                self._keeps_state is not KeepsStateForNextUse.NECESSARY)
 
-        if f_init_not_connected and no_settings_in:
-            do_reuse = self._first_run
-            self._first_run = False
-        else:
-            for message in self._f_init_cache.values():
-                if isinstance(message.data, ClosePort):
-                    do_reuse = False
+        if do_implicit_checkpoint:
+            if self._trigger_manager.should_save_final_snapshot(
+                    do_reuse, self.__f_init_max_timestamp):
+                # store a None instead of a Message
+                self._save_snapshot(None, True, self.__f_init_max_timestamp)
 
         if not do_reuse:
             self.__close_ports()
             self._communicator.shutdown()
             self._deregister()
             self.__manager.close()
+
+        self._api_guard.reuse_instance_done(do_reuse)
         return do_reuse
 
     def error_shutdown(self, message: str) -> None:
@@ -210,7 +272,7 @@ def get_setting(self, name: str, typ: Optional[str] = None
                 self._instance_name(), Reference(name), typ)
 
     def list_ports(self) -> Dict[Operator, List[str]]:
-        """Returns a description of the ports that this CE has.
+        """Returns a description of the ports that this Instance has.
 
         Note that the result has almost the same format as the port
         declarations you pass when making an Instance. The only
@@ -218,9 +280,9 @@ def list_ports(self) -> Dict[Operator, List[str]]:
         even if the port is a vector port.
 
         Returns:
-            A dictionary, indexed by Operator, containing lists of
-            port names. Operators with no associated ports are not
-            included.
+            A dictionary, indexed by :external:py:class:`~ymmsl.Operator`,
+            containing lists of port names. Operators with no associated ports
+            are not included.
         """
         return self._communicator.list_ports()
 
@@ -273,7 +335,8 @@ def get_port_length(self, port: str) -> int:
         Args:
             port: The name of the port to measure.
 
-        Raises: RuntimeError if this is a scalar port.
+        Raises:
+            RuntimeError: If this is a scalar port.
         """
         return self._communicator.get_port(port).get_length()
 
@@ -281,7 +344,7 @@ def set_port_length(self, port: str, length: int) -> None:
         """Resizes the port to the given length.
 
         You should check whether the port is resizable using
-        `is_resizable()` first; whether it is depends on how this
+        :meth:`is_resizable()` first; whether it is depends on how this
         component is wired up, so you should check.
 
         Args:
@@ -298,7 +361,7 @@ def send(self, port_name: str, message: Message,
         """Send a message to the outside world.
 
         Sending is non-blocking, a copy of the message will be made
-        and stored until the receiver is ready to receive it.
+        and stored in memory until the receiver is ready to receive it.
 
         Args:
             port_name: The port on which this message is to be sent.
@@ -310,7 +373,9 @@ def send(self, port_name: str, message: Message,
             message = copy(message)
             message.settings = self._settings_manager.overlay
 
-        self._communicator.send_message(port_name, message, slot)
+        self._communicator.send_message(
+                port_name, message, slot,
+                self._trigger_manager.checkpoints_considered_until())
 
     def receive(self, port_name: str, slot: Optional[int] = None,
                 default: Optional[Message] = None
@@ -381,14 +446,211 @@ def receive_with_settings(
         """
         return self.__receive_message(port_name, slot, default, True)
 
+    def snapshots_enabled(self) -> bool:
+        """Check if the current workflow has snapshots enabled.
+
+        When snapshots are not enabled, all calls to
+        :meth:`should_save_snapshot` and :meth:`should_save_final_snapshot` will
+        return False.
+
+        Returns:
+            True iff checkpoint rules are defined in the workflow yMMSL.
+        """
+        return self._trigger_manager.snapshots_enabled()
+
+    def resuming(self) -> bool:
+        """Check if this instance is resuming from a snapshot.
+
+        Must be used by submodels that implement the checkpointing API. You'll
+        get a RuntimeError when not calling this method in an iteration of the
+        reuse loop.
+
+        This method returns True for the first iteration of the reuse loop after
+        resuming from a previously taken snapshot. When resuming from a
+        snapshot, the submodel must load its state from the snapshot as returned
+        by :meth:`load_snapshot`.
+
+        Returns:
+            True iff the submodel must resume from a snapshot instead of the
+            usual F_INIT step during this iteration of the reuse loop.
+        """
+        self._api_guard.verify_resuming()
+        self._api_guard.resuming_done(self._do_resume)
+        return self._do_resume
+
+    def should_init(self) -> bool:
+        """Check if this instance should initialize.
+
+        Must be used by submodels that implement the checkpointing API.
+
+        When resuming from a previous snapshot, instances need not always
+        execute the F_INIT phase of the submodel execution loop. Use this method
+        before attempting to receive data on F_INIT ports.
+
+        Returns:
+            True if the submodel must execute the F_INIT step, False otherwise.
+        """
+        self._api_guard.verify_should_init()
+        self._api_guard.should_init_done()
+        return self._do_init
+
+    def load_snapshot(self) -> Message:
+        """Load a snapshot.
+
+        Must only be called when :meth:`resuming` returns True.
+
+        Returns:
+            Message object containing the state as saved in a previous run
+            through :meth:`save_snapshot` or :meth:`save_final_snapshot`
+
+        Raises:
+            RuntimeError: if not resuming from a snapshot.
+        """
+        self._api_guard.verify_load_snapshot()
+        result = self._snapshot_manager.load_snapshot()
+        self._api_guard.load_snapshot_done()
+        return result
+
+    def should_save_snapshot(self, timestamp: float) -> bool:
+        """Check if a snapshot should be saved after the S Operator of the
+        submodel.
+
+        This method checks if a snapshot should be saved right now, based on the
+        provided timestamp and passed wallclock time.
+
+        When this method returns True, the submodel must also save a snapshot
+        through :meth:`save_snapshot`. A RuntimeError will be generated when not
+        doing so.
+
+        See also :meth:`should_save_final_snapshot` for the variant that must be
+        called at the end of the reuse loop.
+
+        Args:
+            timestamp: current timestamp of the submodel
+
+        Returns:
+            True iff a snapshot should be taken by the submodel according to the
+            checkpoint rules provided in the ymmsl configuration.
+        """
+        self._api_guard.verify_should_save_snapshot()
+        result = self._trigger_manager.should_save_snapshot(timestamp)
+        self._api_guard.should_save_snapshot_done(result)
+        return result
+
+    def save_snapshot(self, message: Message) -> None:
+        """Save a snapshot after the S Operator of the submodel.
+
+        Before saving a snapshot, you should check using
+        :meth:`should_save_snapshot` if a snapshot should be saved according to
+        the checkpoint rules specified in the ymmsl configuration. You should
+        use the same timestamp in the provided Message object as used to query
+        `should_save_snapshot`.
+
+        Although it is allowed to save a snapshot even when
+        :meth:`should_save_snapshot` returns False, you should avoid this: this
+        situation is not likely to lead to a consistent snapshot over all
+        submodels of the run (and therefore it is not useful to restart from).
+        It could also lead to a lot of snapshot files clogging your file system.
+
+        See also :meth:`save_final_snapshot` for the variant that must be called
+        at the end of the reuse loop.
+
+        Args:
+            message: Message object that is saved as snapshot. The message
+                timestamp attribute should be the same as passed to
+                :meth:`should_save_snapshot`. The data attribute can be used to
+                store the internal state of the submodel.
+        """
+        self._api_guard.verify_save_snapshot()
+        if message is None:
+            raise RuntimeError('Please specify a Message to save as snapshot.')
+        self._save_snapshot(message, False)
+        self._api_guard.save_snapshot_done()
+
+    def should_save_final_snapshot(self, *, apply_overlay: bool = True) -> bool:
+        """Check if a snapshot should be saved at the end of the reuse loop.
+
+        This method checks if a snapshot should be saved now.
+
+        When this method returns True, the submodel must also save a snapshot
+        through :meth:`save_final_snapshot`. A :class:`RuntimeError` will be
+        generated when not doing so.
+
+        See also :meth:`should_save_snapshot` for the variant that may be called
+        inside of a time-integration loop of the submodel.
+
+        .. note::
+            This method will block until it can determine whether a final
+            snapshot should be taken. This means it must also determine if this
+            instance is reused. The optional keyword-only argument
+            `apply_overlay` has the same meaning as for :meth:`reuse_instance`.
+
+        Args:
+            apply_overlay: Whether to apply the received settings
+                overlay or to save it. If you're going to use
+                :meth:`receive_with_settings` on your F_INIT ports, set this to
+                False. If you don't know what that means, just call
+                :meth:`should_save_final_snapshot()` without specifying this and
+                everything will be fine. If it turns out that you did need to
+                specify False, MUSCLE3 will tell you about it in an error
+                message and you can add it still.
+
+        Returns:
+            True iff a final snapshot should be taken by the submodel according
+            to the checkpoint rules provided in the ymmsl configuration.
+        """
+        self._api_guard.verify_should_save_final_snapshot()
+
+        self._do_reuse = self._decide_reuse_instance(apply_overlay)
+        result = self._trigger_manager.should_save_final_snapshot(
+                self._do_reuse, self.__f_init_max_timestamp)
+
+        self._api_guard.should_save_final_snapshot_done(result)
+        return result
+
+    def save_final_snapshot(self, message: Message) -> None:
+        """Save a snapshot at the end of the reuse loop.
+
+        Before saving a snapshot, you should check using
+        :meth:`should_save_final_snapshot` if a snapshot should be saved
+        according to the checkpoint rules specified in the ymmsl configuration.
+
+        Although it is allowed to save a snapshot even when
+        :meth:`should_save_final_snapshot` returns False, you should avoid this:
+        this situation is not likely to lead to a consistent snapshot over all
+        submodels of the run (and therefore it is not useful to restart from).
+        It could also lead to a lot of snapshot files clogging your file system.
+
+        See also :meth:`save_snapshot` for the variant that may be called after
+        each S Operator of the submodel.
+
+        Args:
+            message: Message object that is saved as snapshot. The data
+                attribute can be used to store the internal state of the
+                submodel.
+        """
+        self._api_guard.verify_save_final_snapshot()
+        if message is None:
+            raise RuntimeError('Please specify a Message to save as snapshot.')
+        self._save_snapshot(message, True, self.__f_init_max_timestamp)
+        self._api_guard.save_final_snapshot_done()
+
+    @property
+    def __f_init_max_timestamp(self) -> Optional[float]:
+        """Return max timestamp of pre-received F_INIT messages
+        """
+        return max(
+                (msg.timestamp for msg in self._f_init_cache.values()),
+                default=None)
+
     def _register(self) -> None:
         """Register this instance with the manager.
         """
         register_event = self._profiler.start(ProfileEventType.REGISTER)
         locations = self._communicator.get_locations()
         port_list = self.__list_declared_ports()
-        self.__manager.register_instance(self._instance_name(), locations,
-                                         port_list)
+        self.__manager.register_instance(
+                self._instance_name(), locations, port_list)
         register_event.stop()
         _logger.info('Registered with the manager')
 
@@ -456,6 +718,70 @@ def __set_up_logging(self) -> None:
                                                      self.__manager)
             logging.getLogger().addHandler(self._mmp_handler)
 
+    def _decide_reuse_instance(self, apply_overlay: bool) -> bool:
+        """Decide whether and how to reuse the instance.
+
+        This sets self._first_run, self._do_resume and self._do_init, and
+        returns whether to reuse one more time. This is the real top of
+        the reuse loop, and it gets called by reuse_instance and
+        should_save_final_snapshot.
+        """
+        if self._first_run is None:
+            self._first_run = True
+        elif self._first_run:
+            self._first_run = False
+
+        # resume from intermediate
+        if self._first_run and self._snapshot_manager.resuming_from_intermediate():
+            self._do_resume = True
+            self._do_init = False
+            return True
+
+        f_init_connected = self._have_f_init_connections()
+
+        # resume from final
+        if self._first_run and self._snapshot_manager.resuming_from_final():
+            if f_init_connected:
+                got_f_init_messages = self._pre_receive(apply_overlay)
+                self._do_resume = True
+                self._do_init = True
+                return got_f_init_messages
+            else:
+                self._do_resume = False     # unused
+                self._do_init = False       # unused
+                return False
+
+        # fresh start or resuming from implicit snapshot
+        self._do_resume = False
+
+        # simple straight single run without resuming
+        if not f_init_connected:
+            self._do_init = self._first_run
+            return self._first_run
+
+        # not resuming and f_init connected, run while we get messages
+        got_f_init_messages = self._pre_receive(apply_overlay)
+        self._do_init = got_f_init_messages
+        return got_f_init_messages
+
+    def _save_snapshot(
+            self, message: Optional[Message], final: bool,
+            f_init_max_timestamp: Optional[float] = None) -> None:
+        """Save a snapshot to disk and notify manager.
+
+        Args:
+            message: The data to save
+            final: Whether this is a final snapshot or an intermediate
+                one
+            f_init_max_timestamp: Timestamp for final snapshots
+        """
+        triggers = self._trigger_manager.get_triggers()
+        walltime = self._trigger_manager.elapsed_walltime()
+        timestamp = self._snapshot_manager.save_snapshot(
+                message, final, triggers, walltime,
+                f_init_max_timestamp)
+        self._trigger_manager.update_checkpoints(timestamp)
+
     def __receive_message(
             self, port_name: str, slot: Optional[int],
             default: Optional[Message], with_settings: bool
@@ -475,9 +801,10 @@ def __receive_message(
                 if with_settings and msg.settings is None:
                     err_msg = ('If you use receive_with_settings()'
                                ' on an F_INIT port, then you have to'
-                               ' pass False to reuse_instance(),'
-                               ' otherwise the settings will already'
-                               ' have been applied by MUSCLE.')
+                               ' pass apply_overlay=False to reuse_instance() '
+                               ' and should_save_final_snapshot(),'
+                               ' if applicable, otherwise the settings will'
+                               ' already have been applied by MUSCLE.')
                     self.__shutdown(err_msg)
                     raise RuntimeError(err_msg)
             else:
@@ -500,7 +827,7 @@ def __receive_message(
                     raise RuntimeError(err_msg)
 
         else:
-            msg = self._communicator.receive_message(
+            msg, saved_until = self._communicator.receive_message(
                     port_name, slot, default)
             if port.is_connected() and not port.is_open(slot):
                 err_msg = (('Port {} was closed while trying to'
@@ -512,6 +839,7 @@ def __receive_message(
                 self.__check_compatibility(port_name, msg.settings)
             if not with_settings:
                 msg.settings = None
+            self._trigger_manager.harmonise_wall_time(saved_until)
         return msg
 
     def __make_full_name(self
@@ -583,6 +911,32 @@ def __check_port(self, port_name: str) -> None:
             self.__shutdown(err_msg)
             raise RuntimeError(err_msg)
 
+    def _have_f_init_connections(self) -> bool:
+        """Checks whether we have connected F_INIT ports.
+
+        This includes muscle_settings_in, and any user-defined ports.
+        """
+        ports = self._communicator.list_ports()
+        f_init_connected = any(
+                [self.is_connected(port)
+                 for port in ports.get(Operator.F_INIT, [])])
+        return f_init_connected or self._communicator.settings_in_connected()
+
+    def _pre_receive(self, apply_overlay: bool) -> bool:
+        """Pre-receives on all ports.
+
+        This includes muscle_settings_in and all user-defined ports.
+
+        Returns:
+            True iff no ClosePort messages were received.
+        """
+        all_ports_open = self.__receive_settings()
+        self.__pre_receive_f_init(apply_overlay)
+        for message in self._f_init_cache.values():
+            if isinstance(message.data, ClosePort):
+                all_ports_open = False
+        return all_ports_open
+
     def __receive_settings(self) -> bool:
         """Receives settings on muscle_settings_in.
 
@@ -590,7 +944,7 @@ def __receive_settings(self) -> bool:
             False iff the port is connnected and ClosePort was received.
         """
         default_message = Message(0.0, None, Settings(), Settings())
-        message = self._communicator.receive_message(
+        message, saved_until = self._communicator.receive_message(
                 'muscle_settings_in', None, default_message)
         if isinstance(message.data, ClosePort):
             return False
@@ -607,6 +961,8 @@ def __receive_settings(self) -> bool:
         for key, value in message.data.items():
             settings[key] = value
         self._settings_manager.overlay = settings
+
+        self._trigger_manager.harmonise_wall_time(saved_until)
         return True
 
     def __pre_receive_f_init(self, apply_overlay: bool) -> None:
@@ -616,12 +972,14 @@ def __pre_receive_f_init(self, apply_overlay: bool) -> None:
         in self._f_init_cache.
         """
         def pre_receive(port_name: str, slot: Optional[int]) -> None:
-            msg = self._communicator.receive_message(port_name, slot)
+            msg, saved_until = self._communicator.receive_message(
+                    port_name, slot)
             self._f_init_cache[(port_name, slot)] = msg
             if apply_overlay:
                 self.__apply_overlay(msg)
                 self.__check_compatibility(port_name, msg.settings)
                 msg.settings = None
+            self._trigger_manager.harmonise_wall_time(saved_until)
 
         self._f_init_cache = dict()
         ports = self._communicator.list_ports()
diff --git a/libmuscle/python/libmuscle/manager/manager.py b/libmuscle/python/libmuscle/manager/manager.py
index c8aa76d4..0ec1da3b 100644
--- a/libmuscle/python/libmuscle/manager/manager.py
+++ b/libmuscle/python/libmuscle/manager/manager.py
@@ -11,6 +11,7 @@
 from libmuscle.manager.mmp_server import MMPServer
 from libmuscle.manager.instance_manager import InstanceManager
 from libmuscle.manager.run_dir import RunDir
+from libmuscle.manager.snapshot_registry import SnapshotRegistry
 from libmuscle.manager.topology_store import TopologyStore
 
 
@@ -42,6 +43,14 @@ def __init__(
         self._logger = Logger(log_dir, log_level)
         self._topology_store = TopologyStore(configuration)
         self._instance_registry = InstanceRegistry()
+        if run_dir is not None:
+            snapshot_dir = run_dir.snapshot_dir()
+        else:
+            snapshot_dir = Path.cwd()
+            if self._configuration.checkpoints:
+                _logger.warning('Checkpoints are configured but no run'
+                                ' directory is provided. Snapshots will be'
+                                ' stored in the current working directory.')
 
         if self._run_dir:
             save_ymmsl(
@@ -57,9 +66,16 @@ def __init__(
         except ValueError:
             pass
 
+        # SnapshotRegistry creates a worker thread, must be created after
+        # instance_manager which forks the process
+        self._snapshot_registry = SnapshotRegistry(
+                configuration, snapshot_dir, self._topology_store)
+        self._snapshot_registry.start()
+
         self._server = MMPServer(
-                self._logger, self._configuration.settings,
-                self._instance_registry, self._topology_store)
+                self._logger, self._configuration,
+                self._instance_registry, self._topology_store,
+                self._snapshot_registry, run_dir)
 
         if self._instance_manager:
             self._instance_manager.set_manager_location(
@@ -91,6 +107,8 @@ def stop(self) -> None:
         """Shuts down the manager."""
         # self._server.stop()
         self._server.stop()
+        self._snapshot_registry.shutdown()
+        self._snapshot_registry.join()
         self._logger.close()
 
     def wait(self) -> bool:
diff --git a/libmuscle/python/libmuscle/manager/mmp_server.py b/libmuscle/python/libmuscle/manager/mmp_server.py
index 0f66690a..d609fce1 100644
--- a/libmuscle/python/libmuscle/manager/mmp_server.py
+++ b/libmuscle/python/libmuscle/manager/mmp_server.py
@@ -1,20 +1,25 @@
 import errno
 import logging
-from typing import Any, cast, Generator, List
+import time
+from typing import Any, Dict, cast, List, Optional
 
 import msgpack
-from ymmsl import Conduit, Identifier, Operator, Port, Reference, Settings
+from ymmsl import (
+        Conduit, Identifier, Operator, Port, Reference, PartialConfiguration,
+        Checkpoints)
 
 from libmuscle.logging import LogLevel
 from libmuscle.manager.instance_registry import (
         AlreadyRegistered, InstanceRegistry)
 from libmuscle.manager.logger import Logger
+from libmuscle.manager.run_dir import RunDir
+from libmuscle.manager.snapshot_registry import SnapshotRegistry
 from libmuscle.manager.topology_store import TopologyStore
 from libmuscle.mcp.protocol import RequestType, ResponseType
 from libmuscle.mcp.tcp_transport_server import TcpTransportServer
 from libmuscle.mcp.transport_server import RequestHandler
+from libmuscle.snapshot import SnapshotMetadata
 from libmuscle.timestamp import Timestamp
-from libmuscle.util import generate_indices, instance_indices
 
 
 _logger = logging.getLogger(__name__)
@@ -35,14 +40,26 @@ def encode_conduit(conduit: Conduit) -> List[str]:
     return [str(conduit.sender), str(conduit.receiver)]
 
 
+def encode_checkpoints(checkpoints: Checkpoints) -> Dict[str, Any]:
+    """Convert a Checkpoins to a MsgPack-compatible value."""
+    return {
+        "at_end": checkpoints.at_end,
+        "wallclock_time": [vars(rule) for rule in checkpoints.wallclock_time],
+        "simulation_time": [vars(rule) for rule in checkpoints.simulation_time]
+    }
+
+
 class MMPRequestHandler(RequestHandler):
     """Handles Manager requests."""
     def __init__(
             self,
             logger: Logger,
-            settings: Settings,
+            configuration: PartialConfiguration,
             instance_registry: InstanceRegistry,
-            topology_store: TopologyStore):
+            topology_store: TopologyStore,
+            snapshot_registry: SnapshotRegistry,
+            run_dir: Optional[RunDir]
+            ) -> None:
         """Create an MMPRequestHandler.
 
         Args:
@@ -52,9 +69,12 @@ def __init__(
             topology_store: Keeps track of how to connect things.
         """
         self._logger = logger
-        self._settings = settings
+        self._configuration = configuration
         self._instance_registry = instance_registry
         self._topology_store = topology_store
+        self._snapshot_registry = snapshot_registry
+        self._run_dir = run_dir
+        self._reference_time = time.monotonic()
 
     def handle_request(self, request: bytes) -> bytes:
         """Handles a manager request.
@@ -80,6 +100,10 @@ def handle_request(self, request: bytes) -> bytes:
             response = self._submit_log_message(*req_args)
         elif req_type == RequestType.SUBMIT_PROFILE_EVENTS.value:
             response = self._submit_profile_events(*req_args)
+        elif req_type == RequestType.SUBMIT_SNAPSHOT.value:
+            response = self._submit_snapshot(*req_args)
+        elif req_type == RequestType.GET_CHECKPOINT_INFO.value:
+            response = self._get_checkpoint_info(*req_args)
 
         return cast(bytes, msgpack.packb(response, use_bin_type=True))
 
@@ -100,9 +124,9 @@ def _register_instance(
                 equals ERROR
         """
         port_objs = [decode_port(p) for p in ports]
+        instance = Reference(instance_id)
         try:
-            self._instance_registry.add(
-                Reference(instance_id), locations, port_objs)
+            self._instance_registry.add(instance, locations, port_objs)
 
             _logger.info(f'Registered instance {instance_id}')
             return [ResponseType.SUCCESS.value]
@@ -153,9 +177,10 @@ def _get_peers(self, instance_id: str) -> Any:
 
         # generate instances
         try:
+            peers = self._topology_store.get_peer_instances(instance)
             instance_locations = {
                     str(peer): self._instance_registry.get_locations(peer)
-                    for peer in self._generate_peer_instances(instance)}
+                    for peer in peers}
         except KeyError as e:
             return [
                     ResponseType.PENDING.value,
@@ -202,7 +227,7 @@ def _get_settings(self) -> Any:
         """
         return [
                 ResponseType.SUCCESS.value,
-                self._settings.as_ordered_dict()]
+                self._configuration.settings.as_ordered_dict()]
 
     def _submit_log_message(
             self, instance_id: str, timestamp: float, level: int, text: str
@@ -236,30 +261,51 @@ def _submit_profile_events(self, events: List[List[Any]]) -> Any:
         """
         return [ResponseType.SUCCESS.value]
 
-    def _generate_peer_instances(
-            self, instance: Reference) -> Generator[Reference, None, None]:
-        """Generates the names of all peer instances of an instance.
+    def _submit_snapshot(
+            self, instance_id: str, snapshot: Dict[str, Any]) -> Any:
+        """Handle a submit snapshot request.
+
+        Returns:
+            A list containing the following values on success:
+
+            status (ResponseType): SUCCESS
+        """
+        snapshot_obj = SnapshotMetadata(**snapshot)
+        instance = Reference(instance_id)
+        self._snapshot_registry.register_snapshot(instance, snapshot_obj)
+        return [ResponseType.SUCCESS.value]
+
+    def _get_checkpoint_info(self, instance_id: str) -> Any:
+        """Get checkpoint info for an instance
 
         Args:
-            instance: The instance whose peers to generate.
+            instance: The instance whose checkpoint info to get
+
+        Returns:
+            A list containing the following values on success:
 
-        Yields:
-            All peer instance identifiers.
+            status (ResponseType): SUCCESS
+            wallclock_reference_time (float): Unix timestamp (in UTC) indicating
+                wallclock time of the start of the workflow.
+            checkpoints (dict): Dictionary encoding a ymmsl.Checkpoints object.
+            resume_path (Optional[str]): Checkpoint filename to resume from.
+            snapshot_directory (Optional[str]): Directory to store instance
+                snapshots.
         """
-        component = instance.without_trailing_ints()
-        indices = instance_indices(instance)
-        dims = self._topology_store.kernel_dimensions[component]
-        all_peer_dims = self._topology_store.get_peer_dimensions(component)
-        for peer, peer_dims in all_peer_dims.items():
-            base = peer
-            for i in range(min(len(dims), len(peer_dims))):
-                base += indices[i]
+        instance = Reference(instance_id)
+        resume = None
+        if instance in self._configuration.resume:
+            resume = str(self._configuration.resume[instance])
+
+        snapshot_directory = None
+        if self._run_dir is not None:
+            snapshot_directory = str(self._run_dir.snapshot_dir(instance))
 
-            if dims >= peer_dims:
-                yield base
-            else:
-                for peer_indices in generate_indices(peer_dims[len(dims):]):
-                    yield base + peer_indices
+        return [ResponseType.SUCCESS.value,
+                time.monotonic() - self._reference_time,
+                encode_checkpoints(self._configuration.checkpoints),
+                resume,
+                snapshot_directory]
 
 
 class MMPServer:
@@ -272,9 +318,11 @@ class MMPServer:
     def __init__(
             self,
             logger: Logger,
-            settings: Settings,
+            configuration: PartialConfiguration,
             instance_registry: InstanceRegistry,
-            topology_store: TopologyStore
+            topology_store: TopologyStore,
+            snapshot_registry: SnapshotRegistry,
+            run_dir: Optional[RunDir]
             ) -> None:
         """Create an MMPServer.
 
@@ -285,13 +333,15 @@ def __init__(
 
         Args:
             logger: Logger to send log messages to
-            settings: Settings component to get settings from
+            configuration: Configuration component to get settings, checkpoints
+                and resumes from
             instance_registry: To register instances with and get
                 peer locations from
             topology_store: To get peers and conduits from
         """
         self._handler = MMPRequestHandler(
-                logger, settings, instance_registry, topology_store)
+                logger, configuration, instance_registry, topology_store,
+                snapshot_registry, run_dir)
         try:
             self._server = TcpTransportServer(self._handler, 9000)
         except OSError as e:
diff --git a/libmuscle/python/libmuscle/manager/run_dir.py b/libmuscle/python/libmuscle/manager/run_dir.py
index 8bb4b91e..186d32e8 100644
--- a/libmuscle/python/libmuscle/manager/run_dir.py
+++ b/libmuscle/python/libmuscle/manager/run_dir.py
@@ -1,4 +1,5 @@
 from pathlib import Path
+from typing import Optional
 
 from ymmsl import Reference
 
@@ -20,6 +21,8 @@ class RunDir:
                 <instance_name[i]>.out
                 <instance_name[i]>.err
                 work_dir/
+                snapshots/
+        snapshots/
     """
     def __init__(self, run_dir: Path) -> None:
         """Create a RunDir managing the given directory.
@@ -57,3 +60,20 @@ def instance_dir(self, name: Reference) -> Path:
         make it.
         """
         return self.path / 'instances' / str(name)
+
+    def snapshot_dir(self, name: Optional[Reference] = None) -> Path:
+        """Return the snapshot directory for the workflow or for an instance.
+
+        Args:
+            name: Name of the instance. May be None to get the workflow snapshot
+                directory.
+
+        Returns:
+            The path to the snapshot directory
+        """
+        if name is None:
+            path = self.path / 'snapshots'
+        else:
+            path = self.instance_dir(name) / 'snapshots'
+        path.mkdir(parents=True, exist_ok=True)
+        return path
diff --git a/libmuscle/python/libmuscle/manager/snapshot_registry.py b/libmuscle/python/libmuscle/manager/snapshot_registry.py
new file mode 100644
index 00000000..3883ea9a
--- /dev/null
+++ b/libmuscle/python/libmuscle/manager/snapshot_registry.py
@@ -0,0 +1,588 @@
+from dataclasses import dataclass, field
+from datetime import datetime
+from enum import Flag, auto
+from functools import lru_cache
+from itertools import chain, zip_longest
+from operator import attrgetter
+from pathlib import Path
+from queue import Queue
+from threading import Thread
+from typing import Dict, Optional, Set, FrozenSet, List, Tuple, TypeVar
+
+from ymmsl import (
+        Reference, Model, Identifier, Implementation, save,
+        PartialConfiguration)
+
+from libmuscle.manager.topology_store import TopologyStore
+from libmuscle.snapshot import SnapshotMetadata
+
+
+_MAX_FILE_EXISTS_CHECK = 100
+
+_SnapshotDictType = Dict[Reference, List["SnapshotNode"]]
+_ConnectionType = Tuple[Identifier, Identifier, "_ConnectionInfo"]
+_QueueItemType = Optional[Tuple[Reference, SnapshotMetadata]]
+_T = TypeVar("_T")
+
+# this snapshot is used as a placeholder for restarting from scratch
+_NULL_SNAPSHOT = SnapshotMetadata(["Instance start"], 0, 0, None, {}, True, '')
+
+
+def safe_get(lst: List[_T], index: int, default: _T) -> _T:
+    """Get an item from the list, returning default when it does not exist.
+
+    Args:
+        lst: List to get the item from
+        index: Which item to get, should be >= 0
+        default: Value to return when hitting an IndexError
+    """
+    try:
+        return lst[index]
+    except IndexError:
+        return default
+
+
+class _ConnectionInfo(Flag):
+    SELF_IS_SENDING = auto()
+    SELF_IS_VECTOR = auto()
+    PEER_IS_VECTOR = auto()
+
+
+def calc_consistency(
+        num1: int, num2: int, first_is_sent: bool, num2_is_restart: bool
+        ) -> bool:
+    """Calculate consistency of message counts.
+
+    Args:
+        num1: message count of instance 1
+        num2: message count of instance 2
+        first_is_sent: True iff instance 1 is sending messages over this conduit
+        num2_is_restart: True iff the snapshot of num2 is a full restart
+
+    Returns:
+        True iff the two message counts are consistent
+    """
+    return (num1 == num2 or                             # strong
+            num1 + 1 == num2 and first_is_sent or       # weak (1 = sent)
+            # weak (2 = sent) - only allow if num2 is not a restart
+            num2 + 1 == num1 and not first_is_sent and not num2_is_restart)
+
+
+def calc_consistency_list(
+        num1: List[int], num2: List[int], first_is_sent: bool,
+        num2_is_restart: bool) -> bool:
+    """Calculate consistency of message counts.
+
+    Args:
+        num1: message count of instance 1
+        num2: message count of instance 2
+        first_is_sent: True iff instance 1 is sending messages over this conduit
+        num2_is_restart: True iff the snapshot of num2 is a full restart
+
+    Returns:
+        True iff the two message counts are consistent
+    """
+    if first_is_sent:
+        allow_weak = True
+        slot_iter = zip_longest(num1, num2, fillvalue=0)
+    else:
+        allow_weak = not num2_is_restart
+        slot_iter = zip_longest(num2, num1, fillvalue=0)
+    return all(slot_sent == slot_received or                    # strong
+               slot_sent + 1 == slot_received and allow_weak    # weak
+               for slot_sent, slot_received in slot_iter)
+
+
+@dataclass
+class SnapshotNode:
+    """Represents a node in the snapshot graph.
+
+    Attributes:
+        num: The number of the snapshot. Unique for this instance. Later
+            snapshots always have a higher num.
+        instance: Which instance this is a snapshot of.
+        snapshot: The snapshot metadata reported by the instance.
+        peers: The set of peers that the instance is connected to.
+        consistent_peers: Keeps track of snapshots per peer that are consistent
+            with this one.
+    """
+    num: int
+    instance: Reference
+    snapshot: SnapshotMetadata
+    peers: FrozenSet[Reference]
+    consistent_peers: Dict[Reference, List["SnapshotNode"]] = field(
+            default_factory=dict, repr=False)
+
+    def __hash__(self) -> int:
+        return object.__hash__(self)
+
+    @property
+    def consistent(self) -> bool:
+        """Returns True iff there is a consistent checkpoint with all peers.
+        """
+        return self.consistent_peers.keys() == self.peers
+
+    def do_consistency_check(
+            self,
+            peer_node: "SnapshotNode",
+            connections: List[_ConnectionType]) -> bool:
+        """Check if the snapshot of the peer is consistent with us.
+
+        When the peer snapshot is consistent, adds it to our list of consistent
+        peer snapshots (in :attribute:`consistent_peers`) and vice versa.
+
+        Args:
+            peer_node: Snapshot of one of our peers
+            connections: All connections from our instance to the peer instance
+
+        Returns:
+            True iff the peer snapshot is consistent with ours.
+        """
+        i_snapshot = self.snapshot
+        p_snapshot = peer_node.snapshot
+        peer_is_restart = p_snapshot is _NULL_SNAPSHOT
+        for connection in connections:
+            i_port, p_port, conn = connection
+            is_sending = bool(conn & _ConnectionInfo.SELF_IS_SENDING)
+            i_msg_counts = i_snapshot.port_message_counts.get(str(i_port), [])
+            p_msg_counts = p_snapshot.port_message_counts.get(str(p_port), [])
+            if conn & _ConnectionInfo.SELF_IS_VECTOR:
+                slot = int(peer_node.instance[-1])
+                consistent = calc_consistency(
+                        safe_get(i_msg_counts, slot, 0),
+                        safe_get(p_msg_counts, 0, 0),
+                        is_sending, peer_is_restart)
+            elif conn & _ConnectionInfo.PEER_IS_VECTOR:
+                slot = int(self.instance[-1])
+                consistent = calc_consistency(
+                        safe_get(i_msg_counts, 0, 0),
+                        safe_get(p_msg_counts, slot, 0),
+                        is_sending, peer_is_restart)
+            else:
+                consistent = calc_consistency_list(
+                        i_msg_counts, p_msg_counts, is_sending, peer_is_restart)
+            if not consistent:
+                return False
+        self.consistent_peers.setdefault(
+                peer_node.instance, []).append(peer_node)
+        peer_node.consistent_peers.setdefault(
+                self.instance, []).append(self)
+        return True
+
+
+class SnapshotRegistry(Thread):
+    """Registry of all snapshots taken by instances.
+
+    Current snapshots are stored in a graph. Every node represents a snapshot
+    taken by an instance (see :class:`SnapshotNode`). When snapshots from peer
+    instances are consistent, the nodes are connected to each other.
+
+    This class manages the snapshot nodes. New snapshots are registered through
+    :meth:`register_snapshot`.
+    """
+
+    def __init__(
+            self, config: PartialConfiguration, snapshot_folder: Path,
+            topology_store: TopologyStore) -> None:
+        """Create a snapshot graph using provided configuration.
+
+        Args:
+            config: ymmsl configuration describing the workflow.
+        """
+        super().__init__(name='SnapshotRegistry')
+
+        if config.model is None or not isinstance(config.model, Model):
+            raise ValueError('The yMMSL experiment description does not'
+                             ' contain a (complete) model section, so there'
+                             ' is nothing to run!')
+        self._configuration = config
+        self._model = config.model
+        self._snapshot_folder = snapshot_folder
+        self._topology_store = topology_store
+
+        self._queue = Queue()               # type: Queue[_QueueItemType]
+        self._snapshots = {}                # type: _SnapshotDictType
+
+        self._instances = set()             # type: Set[Reference]
+        for component in config.model.components:
+            self._instances.update(component.instances())
+
+        # Create snapshot nodes for starting from scratch
+        for instance in self._instances:
+            self.register_snapshot(instance, _NULL_SNAPSHOT)
+
+    def register_snapshot(
+            self, instance: Reference, snapshot: SnapshotMetadata) -> None:
+        """Register a new snapshot.
+
+        Args:
+            instance: The instance that created the snapshot
+            snapshot: Metadata describing the snapshot
+        """
+        self._queue.put((instance, snapshot))
+
+    def run(self) -> None:
+        """Code executed in a separate thread
+        """
+        while True:
+            item = self._queue.get()
+            if item is None:
+                return
+            self._add_snapshot(*item)
+
+    def shutdown(self) -> None:
+        """Stop the snapshot registry thread
+        """
+        self._queue.put(None)
+
+    def _add_snapshot(
+            self, instance: Reference, snapshot: SnapshotMetadata) -> None:
+        """Register a new snapshot.
+
+        Args:
+            instance: The instance that created the snapshot
+            snapshot: Metadata describing the snapshot
+        """
+        stateful_peers = self._get_peers(instance)
+
+        i_snapshots = self._snapshots.setdefault(instance, [])
+        # get next number of the snapshot
+        num = 1 if not i_snapshots else i_snapshots[-1].num + 1
+        snapshotnode = SnapshotNode(num, instance, snapshot, stateful_peers)
+        i_snapshots.append(snapshotnode)
+
+        # check consistency with all peers
+        for peer in stateful_peers:
+            for peer_snapshot in self._snapshots.get(peer, []):
+                snapshotnode.do_consistency_check(
+                        peer_snapshot, self._get_connections(instance, peer))
+
+        # finally, check if this snapshotnode is now part of a workflow snapshot
+        if snapshot is not _NULL_SNAPSHOT:
+            self._save_workflow_snapshot(snapshotnode)
+
+    def _save_workflow_snapshot(self, snapshotnode: SnapshotNode) -> None:
+        """Save snapshot if a workflow snapshot exists with the provided node.
+
+        Args:
+            snapshotnode: The snapshot node that must be part of the workflow
+                snapshot.
+        """
+        workflow_snapshots = self._get_workflow_snapshots(snapshotnode)
+        for workflow_snapshot in workflow_snapshots:
+            self._write_snapshot_ymmsl(workflow_snapshot)
+        self._cleanup_snapshots(workflow_snapshots)
+
+    def _get_workflow_snapshots(
+            self, snapshot: SnapshotNode) -> List[List[SnapshotNode]]:
+        """Return all workflow snapshots which contain the provided node.
+
+        Args:
+            snapshotnode: The snapshot node that must be part of the workflow
+                snapshot.
+
+        Returns:
+            List of workflow snapshots. Each workflow snapshot is a list of
+            instance snapshot nodes.
+        """
+        if not snapshot.consistent:
+            return []
+
+        # Instances that don't have a snapshot node chosen yet:
+        instances_to_cover = list(self._instances - {snapshot.instance})
+        # Allowed snapshots per instance. This is updated during the heuristic
+        # to further restrict the sets of snapshots as peer snapshots are
+        # selected.
+        # First restriction is that the snapshots have to be locally consistent.
+        allowed_snapshots = {}  # type: Dict[Reference, FrozenSet[SnapshotNode]]
+        for instance in instances_to_cover:
+            allowed_snapshots[instance] = frozenset(
+                    i_snapshot
+                    for i_snapshot in self._snapshots.get(instance, [])
+                    if i_snapshot.consistent)
+            if not allowed_snapshots[instance]:
+                # there cannot be a workflow snapshot if this instance has no
+                # consistent snapshot nodes
+                return []
+        instance = snapshot.instance
+        allowed_snapshots[instance] = frozenset({snapshot})
+
+        def num_allowed_snapshots(instance: Reference) -> int:
+            """Get number of allowed snapshots at this point for this instance.
+
+            The allowed snapshots are those that are consistent with all
+            selected snapshots at this point in the heuristic.
+            """
+            return len(allowed_snapshots[instance])
+
+        # Do a full, depth-first search for all workflow snapshots
+        # ========================================================
+
+        workflow_snapshots = []
+        selected_snapshots = [snapshot]
+        # This stack stores history of allowed_snapshots and enables roll back
+        stack = []  # type: List[Dict[Reference, FrozenSet[SnapshotNode]]]
+
+        # Update allowed_snapshots for peers of the selected snapshot
+        for peer, snapshots in snapshot.consistent_peers.items():
+            intersection = allowed_snapshots[peer].intersection(snapshots)
+            if not intersection:
+                return []
+            allowed_snapshots[peer] = intersection
+
+        while True:
+            # 1. Select most constrained instance
+            #
+            # Note: we're only interested in the instance with the least allowed
+            # snapshots. Better performance may be possible by not doing a full
+            # sort, but it should be tested. Expectation is that
+            # instances_to_cover remains mostly sorted (as the only counts that
+            # are changing are for peers of the previous selected instance).
+            # Python's sort algorithm is O(N) when the list is already sorted
+            # (which is the same as max()).
+            #
+            # We cannot use a priority queue (heapq) because
+            # num_allowed_snapshots is changing every iteration.
+            instances_to_cover.sort(key=num_allowed_snapshots, reverse=True)
+            instance = instances_to_cover.pop()
+
+            # 2. Select the oldest snapshot of this instance
+            snapshot = min(allowed_snapshots[instance], key=attrgetter('num'))
+            selected_snapshots.append(snapshot)
+            # A shallow copy is ok: the values are immutable frozensets
+            stack.append(allowed_snapshots.copy())
+
+            # 3. Update allowed snapshots based on the newly selected
+            allowed_snapshots[instance] = frozenset({snapshot})
+            for peer, snapshots in snapshot.consistent_peers.items():
+                intersection = allowed_snapshots[peer].intersection(snapshots)
+                if not intersection:
+                    break  # roll back
+                allowed_snapshots[peer] = intersection
+            else:
+                # 4. Selected snapshot is okay to explore further
+                if instances_to_cover:
+                    # 4a. There are still instance to cover, return to the start
+                    #     of the while loop.
+                    continue
+                # 4b. We have found a complete workflow snapshot
+                workflow_snapshots.append(selected_snapshots.copy())
+                # Next: perform a roll-back to continue the search
+
+            # 5. Roll back
+            # stop when selected_snapshots only contains the one we forced to be
+            # part of the workflow snapshot
+            while len(selected_snapshots) > 1:
+                snapshot = selected_snapshots.pop()
+                instance = snapshot.instance
+                instances_to_cover.append(instance)
+                allowed_snapshots = stack.pop()
+                intersection = allowed_snapshots[instance] - {snapshot}
+                allowed_snapshots[instance] = intersection
+                if intersection:
+                    # We have a valid next snapshot to try for this instance
+                    break
+                # No allowed_snapshots, try another roll back
+            else:
+                # Exhausted all roll back possibilities, so we are done now
+                return workflow_snapshots
+
+    def _write_snapshot_ymmsl(
+            self, selected_snapshots: List[SnapshotNode]) -> None:
+        """Write the snapshot ymmsl file to the snapshot folder.
+
+        Args:
+            selected_snapshots: All snapshot nodes of the workflow snapshot.
+        """
+        now = datetime.now()
+        config = self._generate_snapshot_config(selected_snapshots, now)
+        time = now.strftime('%Y%m%d_%H%M%S')
+        for i in range(_MAX_FILE_EXISTS_CHECK):
+            if i:
+                snapshot_filename = f'snapshot_{time}_{i}.ymmsl'
+            else:
+                snapshot_filename = f'snapshot_{time}.ymmsl'
+            savepath = self._snapshot_folder / snapshot_filename
+            if not savepath.exists():
+                save(config, savepath)
+                return
+        raise RuntimeError('Could not find an available filename for storing'
+                           f' the next workflow snapshot: {savepath} already'
+                           ' exists.')
+
+    def _generate_snapshot_config(
+                self, selected_snapshots: List[SnapshotNode], now: datetime
+                ) -> PartialConfiguration:
+        """Generate ymmsl configuration for snapshot file
+        """
+        selected_snapshots.sort(key=attrgetter('instance'))
+        resume = {}
+        for node in selected_snapshots:
+            if node.snapshot is not _NULL_SNAPSHOT:
+                # Only store resume information when it is an actual snapshot
+                # created by the instance. Otherwise the instance can just be
+                # restarted from the beginning.
+                resume[node.instance] = Path(node.snapshot.snapshot_filename)
+        description = self._generate_description(selected_snapshots, now)
+        return PartialConfiguration(resume=resume, description=description)
+
+    def _generate_description(
+            self, selected_snapshots: List[SnapshotNode], now: datetime) -> str:
+        """Generate a human-readable description of the workflow snapshot.
+        """
+        triggers = {}   # type: Dict[str, List[str]]
+        component_info = []
+        max_instance_len = len('Instance ')
+        for node in selected_snapshots:
+            for trigger in node.snapshot.triggers:
+                triggers.setdefault(trigger, []).append(str(node.instance))
+            component_info.append((
+                    str(node.instance),
+                    f'{node.snapshot.timestamp:<11.6g}',
+                    f'{node.snapshot.wallclock_time:<11.6g}',
+                    ("Intermediate", "Final")[node.snapshot.is_final_snapshot]))
+            max_instance_len = max(max_instance_len, len(str(node.instance)))
+        instance_with_padding = 'Instance'.ljust(max_instance_len)
+        component_table = [
+                f'{instance_with_padding} t           Wallclock time  Type',
+                f'{"-" * (max_instance_len + 41)}']
+        component_table += [
+                f'{name.ljust(max_instance_len)} {timestamp} {walltime}'
+                f'     {typ}'
+                for name, timestamp, walltime, typ in component_info]
+        return (f'Workflow snapshot for {self._model.name}'
+                f' taken on {now.strftime("%Y-%m-%d %H:%M:%S")}.\n'
+                'Snapshot triggers:\n' +
+                '\n'.join(f'- {trigger} ({", ".join(triggers[trigger])})'
+                          for trigger in sorted(triggers)) +
+                '\n\n' +
+                '\n'.join(component_table) + '\n')
+
+    def _cleanup_snapshots(
+            self, workflow_snapshots: List[List[SnapshotNode]]) -> None:
+        """Remove all snapshots that are older than the selected snapshots.
+
+        Args:
+            selected_snapshots: All snapshot nodes of a workflow snapshot
+        """
+        if not workflow_snapshots:
+            return
+
+        # Find the newest snapshots per instance
+        newest_snapshots = {snapshot.instance: snapshot
+                            for snapshot in workflow_snapshots[0]}
+        for workflow_snapshot in workflow_snapshots[1:]:
+            for snapshot in workflow_snapshot:
+                if newest_snapshots[snapshot.instance].num < snapshot.num:
+                    newest_snapshots[snapshot.instance] = snapshot
+
+        # Remove all snapshots that are older than the newest snapshots
+        removed_snapshots = set()  # type: Set[SnapshotNode]
+        for snapshot in newest_snapshots.values():
+            all_snapshots = self._snapshots[snapshot.instance]
+            idx = all_snapshots.index(snapshot)
+            self._snapshots[snapshot.instance] = all_snapshots[idx:]
+            removed_snapshots.update(all_snapshots[:idx])
+
+        # Remove all references in SnapshotNode.peer_snapshot to the snapshots
+        # that are cleaned up
+        for snapshot in removed_snapshots:
+            for peer_snapshot in chain.from_iterable(
+                    snapshot.consistent_peers.values()):
+                if peer_snapshot in removed_snapshots:
+                    # snapshot is removed anyway, no need to update references
+                    continue
+                # peer_snapshot is still there, remove reference to us
+                peer_snapshot.consistent_peers[snapshot.instance].remove(
+                        snapshot)
+
+    @lru_cache(maxsize=None)
+    def _get_peers(self, instance: Reference) -> FrozenSet[Reference]:
+        """Return the set of peers for the given instance.
+
+        Note: instance is assumed to contain the full index, not just the
+        component name.
+
+        Args:
+            instance: Instance to get peers of.
+
+        Returns:
+            Frozen set with all peer instances (including their index).
+        """
+        return frozenset(self._topology_store.get_peer_instances(instance))
+
+    @lru_cache(maxsize=None)
+    def _get_connections(self, instance: Reference, peer: Reference
+                         ) -> List[_ConnectionType]:
+        """Get the list of connections between instance and peer.
+
+        Args:
+            instance: Instance reference (including index)
+            peer: Peer reference (including index)
+
+        Returns:
+            A list of tuples describing all conduits between instance and peer:
+                instance_port (Reference): the port of instance that is
+                    connected to
+                peer_port (Reference): the port on the peer instance
+                info (_ConnectionInfo): flag describing the connection. The
+                    instance is sending when
+                    ``info & _ConnectionInfo.SELF_IS_SENDING`` and receiving
+                    otherwise. When the instance port is a vector port and the
+                    peer port is a non-vector port, the flag
+                    ``_ConnectionInfo.SELF_IS_VECTOR`` is set. In the reverse
+                    situation the flag ``_ConnectionInfo.PEER_IS_VECTOR`` is
+                    set. When both ports are vector or non-vector, neither flag
+                    is set.
+        """
+        instance_kernel = instance.without_trailing_ints()
+        peer_kernel = peer.without_trailing_ints()
+
+        connected_ports = []  # type: List[_ConnectionType]
+        for conduit in self._model.conduits:
+            if (conduit.sending_component() == instance_kernel and
+                    conduit.receiving_component() == peer_kernel):
+                conn_type = _ConnectionInfo.SELF_IS_SENDING
+            elif (conduit.receiving_component() == instance_kernel and
+                    conduit.sending_component() == peer_kernel):
+                conn_type = _ConnectionInfo(0)
+            else:
+                continue
+            instance_ndim = (len(instance) - len(instance_kernel))
+            peer_ndim = (len(peer) - len(peer_kernel))
+            if instance_ndim < peer_ndim:
+                conn_type |= _ConnectionInfo.SELF_IS_VECTOR
+            if instance_ndim > peer_ndim:
+                conn_type |= _ConnectionInfo.PEER_IS_VECTOR
+            # we cannot distinguish scalar-scalar vs. vector-vector
+            # but it does not matter for this logic :)
+            if conn_type & _ConnectionInfo.SELF_IS_SENDING:
+                connected_ports.append((
+                        conduit.sending_port(),
+                        conduit.receiving_port(),
+                        conn_type))
+            else:
+                connected_ports.append((
+                        conduit.receiving_port(),
+                        conduit.sending_port(),
+                        conn_type))
+        return connected_ports
+
+    @lru_cache(maxsize=None)
+    def _implementation(self, kernel: Reference) -> Optional[Implementation]:
+        """Return the implementation of a kernel.
+
+        Args:
+            kernel: The kernel to get the implementation for.
+
+        Returns:
+            Implementation for the kernel, or None if not provided in the
+            configuration.
+        """
+        implementation = None
+        for component in self._model.components:
+            if component.name == kernel:
+                implementation = component.implementation
+        if implementation in self._configuration.implementations:
+            return self._configuration.implementations[implementation]
+        return None
diff --git a/libmuscle/python/libmuscle/manager/topology_store.py b/libmuscle/python/libmuscle/manager/topology_store.py
index c8e3f864..f3dd1fd3 100644
--- a/libmuscle/python/libmuscle/manager/topology_store.py
+++ b/libmuscle/python/libmuscle/manager/topology_store.py
@@ -1,4 +1,5 @@
 from typing import Dict, List
+from libmuscle.util import generate_indices, instance_indices
 
 from ymmsl import Conduit, PartialConfiguration, Model, Reference
 
@@ -77,3 +78,30 @@ def get_peer_dimensions(self, kernel_name: Reference
                 snd = conduit.sending_component()
                 ret[snd] = self.kernel_dimensions[snd]
         return ret
+
+    def get_peer_instances(self, instance: Reference) -> List[Reference]:
+        """Generates the names of all peer instances of an instance.
+
+        Args:
+            instance: The instance whose peers to generate.
+
+        Returns:
+            All peer instance identifiers.
+        """
+        component = instance.without_trailing_ints()
+        indices = instance_indices(instance)
+        dims = self.kernel_dimensions[component]
+        all_peer_dims = self.get_peer_dimensions(component)
+
+        peers = []
+        for peer, peer_dims in all_peer_dims.items():
+            base = peer
+            for i in range(min(len(dims), len(peer_dims))):
+                base += indices[i]
+
+            if dims >= peer_dims:
+                peers.append(base)
+            else:
+                for peer_indices in generate_indices(peer_dims[len(dims):]):
+                    peers.append(base + peer_indices)
+        return peers
diff --git a/libmuscle/python/libmuscle/mcp/protocol.py b/libmuscle/python/libmuscle/mcp/protocol.py
index 1e79a11d..5d1217ed 100644
--- a/libmuscle/python/libmuscle/mcp/protocol.py
+++ b/libmuscle/python/libmuscle/mcp/protocol.py
@@ -20,6 +20,8 @@ class RequestType(Enum):
     GET_SETTINGS = 4
     SUBMIT_LOG_MESSAGE = 5
     SUBMIT_PROFILE_EVENTS = 6
+    SUBMIT_SNAPSHOT = 7
+    GET_CHECKPOINT_INFO = 8
 
     # MUSCLE Peer Protocol
     GET_NEXT_MESSAGE = 21
diff --git a/libmuscle/python/libmuscle/mmp_client.py b/libmuscle/python/libmuscle/mmp_client.py
index fd236b52..eed4d99a 100644
--- a/libmuscle/python/libmuscle/mmp_client.py
+++ b/libmuscle/python/libmuscle/mmp_client.py
@@ -1,14 +1,19 @@
+import dataclasses
+from pathlib import Path
 from random import uniform
 from time import perf_counter, sleep
-from typing import Any, Dict, Iterable, List, Tuple
+from typing import Any, Dict, Iterable, List, Optional, Tuple
 
 import msgpack
-from ymmsl import Conduit, Operator, Port, Reference, Settings
+from ymmsl import (
+        Conduit, Operator, Port, Reference, Settings, Checkpoints,
+        CheckpointRule, CheckpointRangeRule, CheckpointAtRule)
 
 from libmuscle.mcp.protocol import RequestType, ResponseType
 from libmuscle.mcp.tcp_transport_client import TcpTransportClient
 from libmuscle.profiling import ProfileEvent
 from libmuscle.logging import LogMessage
+from libmuscle.snapshot import SnapshotMetadata
 
 
 CONNECTION_TIMEOUT = 300
@@ -16,6 +21,9 @@
 PEER_INTERVAL_MIN = 5.0
 PEER_INTERVAL_MAX = 10.0
 
+_CheckpointInfoType = Tuple[
+        float, Checkpoints, Optional[Path], Optional[Path]]
+
 
 def encode_operator(op: Operator) -> str:
     """Convert an Operator to a MsgPack-compatible value."""
@@ -45,6 +53,46 @@ def encode_profile_event(event: ProfileEvent) -> Any:
             event.message_size]
 
 
+def decode_checkpoint_rule(rule: Dict[str, Any]) -> CheckpointRule:
+    """Decode a checkpoint rule from a MsgPack-compatible value."""
+    if rule.keys() == {'at'}:
+        return CheckpointAtRule(**rule)
+    if rule.keys() == {'start', 'stop', 'every'}:
+        return CheckpointRangeRule(**rule)
+    raise ValueError(f'Cannot convert {rule} to a checkpoint rule.')
+
+
+def decode_checkpoint_info(
+        elapsed_time: float,
+        checkpoints_dict: Dict[str, Any],
+        resume: Optional[str],
+        snapshot_dir: Optional[str]
+        ) -> _CheckpointInfoType:
+    """Decode checkpoint info from a MsgPack-compatible value.
+
+    Args:
+        elapsed_time: current elapsed time according to the manager
+        checkpoints_dict: checkpoint definitions from the MsgPack
+        resume: path to the snapshot we should resume from, if any
+        snapshot_dir: path to the directory to store new snapshots in
+
+    Returns:
+        elapsed_time: current elapsed time according to the manager
+        checkpoints: checkpoint configuration
+        resume: path to the snapshot we should resume from, if any
+        snapshot_dir: path to the directory to store new snapshots in
+    """
+    checkpoints = Checkpoints(
+            at_end=checkpoints_dict["at_end"],
+            wallclock_time=[decode_checkpoint_rule(rule)
+                            for rule in checkpoints_dict["wallclock_time"]],
+            simulation_time=[decode_checkpoint_rule(rule)
+                             for rule in checkpoints_dict["simulation_time"]])
+    resume_path = None if resume is None else Path(resume)
+    snapshot_path = None if snapshot_dir is None else Path(snapshot_dir)
+    return (elapsed_time, checkpoints, resume_path, snapshot_path)
+
+
 class MMPClient():
     """The client for the MUSCLE Manager Protocol.
 
@@ -93,6 +141,21 @@ def submit_profile_events(self, events: Iterable[ProfileEvent]) -> None:
                 [encode_profile_event(e) for e in events]]
         self._call_manager(request)
 
+    def submit_snapshot_metadata(
+                self, name: Reference, snapshot_metadata: SnapshotMetadata
+                ) -> None:
+        """Send snapshot metadata to the manager.
+
+        Args:
+            name: Name of the instance in the simulation.
+            snapshot_metadata: Snapshot metadata to supply to the manager.
+        """
+        request = [
+                RequestType.SUBMIT_SNAPSHOT.value,
+                str(name),
+                dataclasses.asdict(snapshot_metadata)]
+        self._call_manager(request)
+
     def get_settings(self) -> Settings:
         """Get the central settings from the manager.
 
@@ -103,6 +166,19 @@ def get_settings(self) -> Settings:
         response = self._call_manager(request)
         return Settings(response[1])
 
+    def get_checkpoint_info(self, name: Reference) -> _CheckpointInfoType:
+        """Get the checkpoint info from the manager.
+
+        Returns:
+            elapsed_time: current elapsed time
+            checkpoints: checkpoint configuration
+            resume: path to the resume snapshot
+            snapshot_directory: path to store snapshots
+        """
+        request = [RequestType.GET_CHECKPOINT_INFO.value, str(name)]
+        response = self._call_manager(request)
+        return decode_checkpoint_info(*response[1:])
+
     def register_instance(self, name: Reference, locations: List[str],
                           ports: List[Port]) -> None:
         """Register a component instance with the manager.
@@ -118,7 +194,7 @@ def register_instance(self, name: Reference, locations: List[str],
                 str(name), locations,
                 [encode_port(p) for p in ports]]
         response = self._call_manager(request)
-        if len(response) > 1:
+        if response[0] == ResponseType.ERROR.value:
             raise RuntimeError(
                     f'Error registering instance: {response[1]}')
 
diff --git a/libmuscle/python/libmuscle/mpp_message.py b/libmuscle/python/libmuscle/mpp_message.py
index 8b84aeab..976f9f5e 100644
--- a/libmuscle/python/libmuscle/mpp_message.py
+++ b/libmuscle/python/libmuscle/mpp_message.py
@@ -93,7 +93,7 @@ def _decode_grid(code: int, data: bytes) -> Grid:
             ExtTypeId.GRID_INT64: np.int64,
             ExtTypeId.GRID_FLOAT32: np.float32,
             ExtTypeId.GRID_FLOAT64: np.float64,
-            ExtTypeId.GRID_BOOL: np.bool8}
+            ExtTypeId.GRID_BOOL: np.bool_}
 
     order_map = {
             'fa': 'F',
@@ -151,7 +151,8 @@ class MPPMessage:
     def __init__(self, sender: Reference, receiver: Reference,
                  port_length: Optional[int],
                  timestamp: float, next_timestamp: Optional[float],
-                 settings_overlay: Settings, data: Any
+                 settings_overlay: Settings, message_number: int,
+                 saved_until: float, data: Any
                  ) -> None:
         """Create an MPPMessage.
 
@@ -169,6 +170,9 @@ def __init__(self, sender: Reference, receiver: Reference,
             receiver: The receiving endpoint.
             port_length: Length of the slot, where applicable.
             settings_overlay: The serialised overlay settings.
+            message_number: Sequence number on this conduit.
+            saved_until: Elapsed time until which the sender has
+                processed checkpoints.
             data: The serialised contents of the message.
         """
         # make sure timestamp and next_timestamp are floats
@@ -182,6 +186,8 @@ def __init__(self, sender: Reference, receiver: Reference,
         self.timestamp = timestamp
         self.next_timestamp = next_timestamp
         self.settings_overlay = settings_overlay
+        self.message_number = message_number
+        self.saved_until = saved_until
         if isinstance(data, np.ndarray):
             self.data = Grid(data)
         else:
@@ -202,11 +208,13 @@ def from_bytes(message: bytes) -> 'MPPMessage':
         timestamp = message_dict["timestamp"]
         next_timestamp = message_dict["next_timestamp"]
         settings_overlay = message_dict["settings_overlay"]
+        message_number = message_dict["message_number"]
+        saved_until = message_dict["saved_until"]
 
         data = message_dict["data"]
         return MPPMessage(
                 sender, receiver, port_length, timestamp, next_timestamp,
-                settings_overlay, data)
+                settings_overlay, message_number, saved_until, data)
 
     def encoded(self) -> bytes:
         """Encode the message and return as a bytes buffer.
@@ -218,6 +226,8 @@ def encoded(self) -> bytes:
                 'timestamp': self.timestamp,
                 'next_timestamp': self.next_timestamp,
                 'settings_overlay': self.settings_overlay,
+                'message_number': self.message_number,
+                'saved_until': self.saved_until,
                 'data': self.data
                 }
 
diff --git a/libmuscle/python/libmuscle/port.py b/libmuscle/python/libmuscle/port.py
index beb49249..a6f955a5 100644
--- a/libmuscle/python/libmuscle/port.py
+++ b/libmuscle/python/libmuscle/port.py
@@ -1,8 +1,20 @@
-from typing import List, Optional
+from typing import List, Optional, TypeVar
+
 from ymmsl import Identifier, Operator
 import ymmsl
 
 
+_T = TypeVar("_T")
+
+
+def _extend_list_to_size(lst: List[_T], size: int, padding: _T) -> None:
+    """When lst is smaller than size, extend to size using padding as values
+    """
+    num_extend = size - len(lst)
+    if num_extend > 0:
+        lst += [padding] * num_extend
+
+
 class Port(ymmsl.Port):
     """Represents a gateway to the outside world.
 
@@ -10,10 +22,15 @@ class Port(ymmsl.Port):
     an operator, as well as a set of dimensions that determines the
     valid slot indices for sending or receiving on this port.
 
+    Ports keep track of the amount of messages sent or received on the port.
+    However, the actual incrementing and validation is done in
+    :class:`Communicator`.
+
     Attributes:
         name (Identifier): Name of this port.
         operator (Operator): Operator associated with this port.
     """
+
     def __init__(self, name: str, operator: Operator, is_vector: bool,
                  is_connected: bool, our_ndims: int, peer_dims: List[int]
                  ) -> None:
@@ -68,6 +85,8 @@ def __init__(self, name: str, operator: Operator, is_vector: bool,
             self._is_open = [True]
 
         self._is_resizable = is_vector and (our_ndims == len(peer_dims))
+        self._num_messages = [0] * (self._length or 1)
+        self._is_resuming = [False] * (self._length or 1)
 
     # Note: I'm not sure how this will develop exactly, so this class has some
     # accessors even if those are un-Pythonic; in the future a simple variable
@@ -129,6 +148,11 @@ def set_length(self, length: int) -> None:
         if length != self._length:
             self._length = length
             self._is_open = [True] * self._length
+            # Using extend here to not discard any information about message
+            # numbers between resizes. Note that _num_messages and _is_resuming
+            # may be longer than self._length!
+            _extend_list_to_size(self._num_messages, self._length, 0)
+            _extend_list_to_size(self._is_resuming, self._length, False)
 
     def set_closed(self, slot: Optional[int] = None) -> None:
         """Marks this port as closed.
@@ -137,3 +161,52 @@ def set_closed(self, slot: Optional[int] = None) -> None:
             self._is_open[slot] = False
         else:
             self._is_open = [False]
+
+    def restore_message_counts(self, num_messages: List[int]) -> None:
+        """Restore message counts from a snapshot
+        """
+        self._num_messages = num_messages
+        self._is_resuming = [True] * len(self._num_messages)
+        _extend_list_to_size(self._num_messages, self._length or 1, 0)
+        _extend_list_to_size(self._is_resuming, self._length or 1, False)
+
+    def get_message_counts(self) -> List[int]:
+        """Get a list of message counts for all slots in this port
+        """
+        return self._num_messages.copy()
+
+    def increment_num_messages(self, slot: Optional[int] = None) -> None:
+        """Increment amount of messages sent or received.
+
+        Args:
+            slot: The slot that is sent/received on
+        """
+        self._num_messages[slot or 0] += 1
+        self.set_resumed(slot)
+
+    def get_num_messages(self, slot: Optional[int] = None) -> int:
+        """Get the amount of messages sent or received.
+
+        Args:
+            slot: The slot that is sent/received on
+        """
+        return self._num_messages[slot or 0]
+
+    def is_resuming(self, slot: Optional[int] = None) -> bool:
+        """True when this port has resumed.
+
+        After resumption, each port/slot may discard exactly one message.
+        is_resuming keeps track of this state.
+
+        Args:
+            slot: The slot that is sent/received on
+        """
+        return self._is_resuming[slot or 0]
+
+    def set_resumed(self, slot: Optional[int] = None) -> None:
+        """Mark that this port has resumed and may no longer discard messages.
+
+        Args:
+            slot: The slot that is sent/received on
+        """
+        self._is_resuming[slot or 0] = False
diff --git a/libmuscle/python/libmuscle/snapshot.py b/libmuscle/python/libmuscle/snapshot.py
new file mode 100644
index 00000000..2f86a220
--- /dev/null
+++ b/libmuscle/python/libmuscle/snapshot.py
@@ -0,0 +1,128 @@
+from abc import ABC, abstractmethod
+from dataclasses import dataclass
+from typing import Dict, List, Optional, cast
+
+import msgpack
+from ymmsl import Reference, Settings
+
+from libmuscle.mpp_message import MPPMessage
+from libmuscle import communicator
+
+
+class Snapshot(ABC):
+    """Snapshot data structure.
+
+    This is an abstract base class, implementations are provided by subclasses.
+    """
+    SNAPSHOT_VERSION_BYTE = b'\0'
+
+    def __init__(self,
+                 triggers: List[str],
+                 wallclock_time: float,
+                 port_message_counts: Dict[str, List[int]],
+                 is_final_snapshot: bool,
+                 message: Optional['communicator.Message']) -> None:
+        self.triggers = triggers
+        self.wallclock_time = wallclock_time
+        self.port_message_counts = port_message_counts
+        self.is_final_snapshot = is_final_snapshot
+        self.message = message
+
+    @classmethod
+    @abstractmethod
+    def from_bytes(cls, data: bytes) -> 'Snapshot':
+        """Create a snapshot object from binary data.
+
+        Args:
+            data: binary data representing the snapshot. Note that this must
+                **exclude** the versioning byte.
+        """
+        ...
+
+    @abstractmethod
+    def to_bytes(self) -> bytes:
+        """Convert the snapshot object to binary data.
+
+        Returns:
+            Binary data representing the snapshot. Note that this must
+                **exclude** the versioning byte.
+        """
+        ...
+
+
+class MsgPackSnapshot(Snapshot):
+    """Snapshot stored in messagepack format
+    """
+    SNAPSHOT_VERSION_BYTE = b'1'
+
+    @classmethod
+    def from_bytes(cls, data: bytes) -> 'Snapshot':
+        dct = msgpack.loads(data)
+        return cls(dct['triggers'],
+                   dct['wallclock_time'],
+                   dct['port_message_counts'],
+                   dct['is_final_snapshot'],
+                   cls.bytes_to_message(dct['message']))
+
+    def to_bytes(self) -> bytes:
+        return cast(bytes, msgpack.dumps({
+            'triggers': self.triggers,
+            'wallclock_time': self.wallclock_time,
+            'port_message_counts': self.port_message_counts,
+            'is_final_snapshot': self.is_final_snapshot,
+            'message': self.message_to_bytes(self.message)
+        }))
+
+    @staticmethod
+    def message_to_bytes(message: Optional['communicator.Message']) -> bytes:
+        """Use MPPMessage serializer for serializing the message object
+        """
+        if message is None:
+            return b''
+        settings = Settings()
+        if message.settings is not None:
+            settings = message.settings
+        return MPPMessage(Reference('_'), Reference('_'), None,
+                          message.timestamp, message.next_timestamp,
+                          settings, 0, -1.0, message.data).encoded()
+
+    @staticmethod
+    def bytes_to_message(data: bytes) -> Optional['communicator.Message']:
+        """Use MPPMessage deserializer for serializing the message object
+        """
+        if not data:
+            return None
+        mpp_message = MPPMessage.from_bytes(data)
+        return communicator.Message(mpp_message.timestamp,
+                                    mpp_message.next_timestamp,
+                                    mpp_message.data,
+                                    mpp_message.settings_overlay)
+
+
+@dataclass
+class SnapshotMetadata:
+    """Metadata of a snapshot for sending to the muscle_manager.
+    """
+    triggers: List[str]
+    wallclock_time: float
+    timestamp: float
+    next_timestamp: Optional[float]
+    port_message_counts: Dict[str, List[int]]
+    is_final_snapshot: bool
+    # storing as str, because Path cannot be serialized by msgpack
+    snapshot_filename: str
+
+    @staticmethod
+    def from_snapshot(snapshot: Snapshot, snapshot_filename: str
+                      ) -> 'SnapshotMetadata':
+        """Create snapshot metadata from the given snapshot and filename
+        """
+        return SnapshotMetadata(
+            snapshot.triggers,
+            snapshot.wallclock_time,
+            snapshot.message.timestamp if snapshot.message else float('NaN'),
+            snapshot.message.next_timestamp if snapshot.message else None,
+            snapshot.port_message_counts,
+            snapshot.is_final_snapshot,
+            snapshot_filename
+        )
diff --git a/libmuscle/python/libmuscle/snapshot_manager.py b/libmuscle/python/libmuscle/snapshot_manager.py
new file mode 100644
index 00000000..f756d05d
--- /dev/null
+++ b/libmuscle/python/libmuscle/snapshot_manager.py
@@ -0,0 +1,211 @@
+import logging
+from pathlib import Path
+from typing import cast, List, Optional
+
+from ymmsl import Reference, Operator
+
+from libmuscle.communicator import Communicator, Message
+from libmuscle.mmp_client import MMPClient
+from libmuscle.snapshot import MsgPackSnapshot, Snapshot, SnapshotMetadata
+
+_logger = logging.getLogger(__name__)
+
+_MAX_FILE_EXISTS_CHECK = 10000
+
+# error text for save_snapshot when msg = None
+_NO_MESSAGE_PROVIDED = (
+        'Invalid message provided to `{}`. Please create a Message object to'
+        ' store the state of the instance in a snapshot.')
+
+
+class SnapshotManager:
+    """Manages information on snapshots for the Instance
+
+    Implements the public checkpointing API with handoffs to
+    :class:`TriggerManager` for checkpoint triggers.
+    """
+
+    def __init__(self,
+                 instance_id: Reference,
+                 manager: MMPClient,
+                 communicator: Communicator) -> None:
+        """Create a new snapshot manager
+
+        Args:
+            instance_id: The id of this instance.
+            manager: The client used to submit data to the manager.
+            communicator: The communicator belonging to this instance.
+        """
+        self._instance_id = instance_id
+        # replace identifier[i] by identifier-i to use in snapshot file name
+        # using a dash (-) because that is not allowed in Identifiers
+        self._safe_id = str(instance_id).replace("[", "-").replace("]", "")
+        self._communicator = communicator
+        self._manager = manager
+
+        self._resume_from_snapshot = None   # type: Optional[Snapshot]
+        self._next_snapshot_num = 1
+
+    def prepare_resume(
+            self, resume_snapshot: Optional[Path],
+            snapshot_directory: Optional[Path]) -> Optional[float]:
+        """Apply checkpoint info received from the manager.
+
+        If there is a snapshot to resume from, this loads it and does
+        any resume work that libmuscle should do, including restoring
+        message counts and storing the resumed-from snapshot again as
+        our first snapshot.
+
+        Args:
+            resume_snapshot: Snapshot to resume from (or None if not
+                resuming)
+            snapshot_directory: directory to save snapshots in
+
+        Returns:
+            Time at which the initial snapshot was saved, if resuming.
+        """
+        result = None       # type: Optional[float]
+        self._snapshot_directory = snapshot_directory or Path.cwd()
+        if resume_snapshot is not None:
+            snapshot = self.load_snapshot_from_file(resume_snapshot)
+
+            if snapshot.message is not None:
+                # snapshot.message is None for implicit snapshots
+                self._resume_from_snapshot = snapshot
+                result = snapshot.message.timestamp
+
+            self._communicator.restore_message_counts(
+                snapshot.port_message_counts)
+            # Store a copy of the snapshot in the current run directory
+            path = self.__store_snapshot(snapshot)
+            metadata = SnapshotMetadata.from_snapshot(snapshot, str(path))
+            self._manager.submit_snapshot_metadata(self._instance_id, metadata)
+
+        return result
+
+    def resuming_from_intermediate(self) -> bool:
+        """Check whether we have an intermediate snapshot.
+
+        Doesn't say whether we should resume now, just that we were
+        given an intermediate snapshot to resume from by the manager.
+        """
+        return (
+                self._resume_from_snapshot is not None and
+                not self._resume_from_snapshot.is_final_snapshot)
+
+    def resuming_from_final(self) -> bool:
+        """Check whether we have a final snapshot.
+
+        Doesn't say whether we should resume now, just that we were
+        given an intermediate snapshot to resume from by the manager.
+        """
+        return (
+                self._resume_from_snapshot is not None and
+                self._resume_from_snapshot.is_final_snapshot)
+
+    def load_snapshot(self) -> Message:
+        """Get the Message to resume from.
+        """
+        snapshot = cast(Snapshot, self._resume_from_snapshot)
+        return cast(Message, snapshot.message)
+
+    def save_snapshot(
+            self, msg: Optional[Message], final: bool,
+            triggers: List[str], wallclock_time: float,
+            f_init_max_timestamp: Optional[float] = None,
+            ) -> float:
+        """Save a (final) snapshot.
+
+        Args:
+            msg: Message object representing the snapshot.
+            final: True iff called from save_final_snapshot.
+            triggers: Description of checkpoints that triggered this.
+            wallclock_time: Wallclock time when saving.
+            f_init_max_timestamp: Timestamp for final snapshots.
+
+        Returns:
+            Simulation time at which the snapshot was made
+        """
+        port_message_counts = self._communicator.get_message_counts()
+        if final:
+            # Decrease F_INIT port counts by one: F_INIT messages are already
+            # pre-received, but not yet processed by the user code. Therefore,
+            # the snapshot state should treat these as not-received.
+            all_ports = self._communicator.list_ports()
+            ports = all_ports.get(Operator.F_INIT, [])
+            if self._communicator.settings_in_connected():
+                ports.append('muscle_settings_in')
+            for port_name in ports:
+                new_counts = [i - 1 for i in port_message_counts[port_name]]
+                port_message_counts[port_name] = new_counts
+
+        snapshot = MsgPackSnapshot(
+            triggers, wallclock_time, port_message_counts, final, msg)
+
+        path = self.__store_snapshot(snapshot)
+        metadata = SnapshotMetadata.from_snapshot(snapshot, str(path))
+        self._manager.submit_snapshot_metadata(self._instance_id, metadata)
+
+        timestamp = msg.timestamp if msg is not None else float('-inf')
+        if final and f_init_max_timestamp is not None:
+            # For final snapshots f_init_max_snapshot is the reference time (see
+            # should_save_final_snapshot).
+            timestamp = f_init_max_timestamp
+        return timestamp
+
+    @staticmethod
+    def load_snapshot_from_file(snapshot_location: Path) -> Snapshot:
+        """Load a previously stored snapshot from the filesystem
+
+        Args:
+            snapshot_location: path where the snapshot is stored
+        """
+        _logger.debug(f'Loading snapshot from {snapshot_location}')
+        if not snapshot_location.is_file():
+            raise RuntimeError(f'Unable to load snapshot: {snapshot_location}'
+                               ' is not a file. Please ensure this path exists'
+                               ' and can be read.')
+
+        # TODO: encapsulate I/O errors?
+        with snapshot_location.open('rb') as snapshot_file:
+            version = snapshot_file.read(1)
+            data = snapshot_file.read()
+
+            if version == MsgPackSnapshot.SNAPSHOT_VERSION_BYTE:
+                return MsgPackSnapshot.from_bytes(data)
+            raise RuntimeError('Unable to load snapshot from'
+                               f' {snapshot_location}: unknown version of'
+                               ' snapshot file. Was the file saved with a'
+                               ' different version of libmuscle or'
+                               ' edited?')
+
+    def __store_snapshot(self, snapshot: Snapshot) -> Path:
+        """Store a snapshot on the filesystem
+
+        Args:
+            snapshot: snapshot to store
+
+        Returns:
+            Path where the snapshot is stored
+        """
+        _logger.debug(f'Saving snapshot to {self._snapshot_directory}')
+        for _ in range(_MAX_FILE_EXISTS_CHECK):
+            # Expectation is that muscle_snapshot_directory is empty initially
+            # and we succeed in the first loop. Still wrapping in a for-loop
+            # such that an existing filename doesn't immediately raise an error
+            fname = f"{self._safe_id}_{self._next_snapshot_num}.pack"
+            fpath = self._snapshot_directory / fname
+            self._next_snapshot_num += 1
+            if not fpath.exists():
+                break
+        else:
+            raise RuntimeError('Could not find an available filename for'
+                               f' storing the next snapshot: {fpath} already'
+                               ' exists.')
+        # Opening with mode 'x' since a file with the same name may be created
+        # in the small window between checking above and opening here. It is
+        # better to fail with an error than to overwrite an existing file.
+        with fpath.open('xb') as snapshot_file:
+            snapshot_file.write(snapshot.SNAPSHOT_VERSION_BYTE)
+            snapshot_file.write(snapshot.to_bytes())
+        return fpath
diff --git a/libmuscle/python/libmuscle/util.py b/libmuscle/python/libmuscle/util.py
index 66217c21..f34d8cd9 100644
--- a/libmuscle/python/libmuscle/util.py
+++ b/libmuscle/python/libmuscle/util.py
@@ -1,3 +1,4 @@
+import itertools
 from pathlib import Path
 import sys
 from typing import Generator, List, Optional, cast
@@ -47,32 +48,8 @@ def generate_indices(dims: List[int]) -> Generator[List[int], None, None]:
     Yields:
         Lists of indices, one for each point in the block.
     """
-    index = [0] * len(dims)
-    done = False
-    while not done:
-        yield index
-        done = increment_index(index, dims)
-
-
-def increment_index(index: List[int], dims: List[int]) -> bool:
-    """Increments an index.
-
-    Args:
-        index: The index to be incremented.
-        dims: The dimensions of the block this index is in.
-
-    Returns:
-        True iff the index overflowed and is now all zeros again.
-    """
-    cur = len(index) - 1
-    index[cur] += 1
-    while index[cur] == dims[cur]:
-        index[cur] = 0
-        if cur == 0:
-            return True
-        cur -= 1
-        index[cur] += 1
-    return False
+    for index in itertools.product(*map(range, dims)):
+        yield list(index)
 
 
 def extract_log_file_location(filename: str) -> Optional[Path]:
diff --git a/muscle3/muscle3.py b/muscle3/muscle3.py
index 6682215f..f715d5fd 100644
--- a/muscle3/muscle3.py
+++ b/muscle3/muscle3.py
@@ -1,12 +1,16 @@
 import sys
+from collections import OrderedDict
+from pathlib import Path
 from typing import Sequence
 
 import click
 import ymmsl
-from ymmsl import Identifier, PartialConfiguration
+from ymmsl import PartialConfiguration
 
 
-from libmuscle.planner.planner import Planner, Resources
+from libmuscle.planner.planner import (
+        Planner, Resources, InsufficientResourcesAvailable)
+from libmuscle.snapshot_manager import SnapshotManager
 
 
 _RESOURCES_INCOMPLETE_MODEL = """
@@ -17,13 +21,11 @@
 
 
 @click.group()
-def muscle3():
+def muscle3() -> None:
     """MUSCLE3 command line interface
 
-    In the future, this command will provide various functions for
-    running coupled simulations using MUSCLE3. For now, it does only
-    one thing, which is to calculate the number of cluster nodes
-    needed for a given simulation to run without oversubscribing.
+    This command provides various functions for running coupled simulations
+    using MUSCLE3.
 
     Use muscle3 <command> --help for help with individual commands.
     """
@@ -107,6 +109,51 @@ def resources(
     sys.exit(0)
 
 
+@muscle3.command(short_help='Display details of a stored snapshot')
+@click.argument(
+        'snapshot_files', nargs=-1, required=True, type=click.Path(
+            exists=True, file_okay=True, dir_okay=False, readable=True,
+            allow_dash=True, resolve_path=True, path_type=Path))
+@click.option(
+        '-d', '--data', is_flag=True,
+        help='Display stored data. Note this may result in a lot of output!')
+@click.option(
+        '-v', '--verbose', is_flag=True, help='Display more metadata.')
+def snapshot(
+        snapshot_files: Sequence[Path], data: bool, verbose: bool) -> None:
+    """Display information about stored snapshots.
+
+    Per provided snapshot, display metadata. Stored data can also be output by
+    supplying the '-d' or '--data' flags. Note that this may result in a lot of
+    data displayed.
+    """
+    for file in snapshot_files:
+        snapshot = SnapshotManager.load_snapshot_from_file(file)
+        click.echo(f'Snapshot at {file}:')
+        typ = 'Final' if snapshot.is_final_snapshot else 'Intermediate'
+        properties = OrderedDict([
+            ('Snapshot type', typ),
+            ('Snapshot timestamp',
+             snapshot.message.timestamp if snapshot.message else float('-inf')),
+            ('Snapshot wallclock time', snapshot.wallclock_time),
+            ('Snapshot triggers', snapshot.triggers),
+        ])
+        if verbose:
+            properties.update([
+                ('Internal: Port message counts', snapshot.port_message_counts),
+            ])
+        for prop_name, prop_value in properties.items():
+            click.secho(f'{prop_name}: ', nl=False, bold=True)
+            click.echo(prop_value)
+        if data:
+            click.secho('Snapshot data:', bold=True)
+            if snapshot.message is not None:
+                click.echo(snapshot.message.data)
+            else:
+                click.secho("No data available", italic=True)
+        click.echo()
+
+
 def _load_ymmsl_files(ymmsl_files: Sequence[str]) -> PartialConfiguration:
     """Loads and merges yMMSL files."""
     configuration = PartialConfiguration()
diff --git a/tox.ini b/tox.ini
index 6f71615f..020ee9e2 100644
--- a/tox.ini
+++ b/tox.ini
@@ -8,7 +8,7 @@ deps =
     flake8<6
     pytest
     pytest-cov
-    git+https://github.com/multiscale/ymmsl-python.git@feature/multicast#egg=ymmsl
+    git+https://github.com/multiscale/ymmsl-python.git@develop#egg=ymmsl
 
 passenv =
     MUSCLE_TEST_PYTHON_ONLY
@@ -17,7 +17,7 @@ passenv =
 
 commands =
     mypy
-    pytest
+    pytest {posargs}
     flake8 libmuscle/python/libmuscle integration_test
 
 [gh-actions]


Test patch:
diff --git a/integration_test/conftest.py b/integration_test/conftest.py
index ad59842a..a6c70b1e 100644
--- a/integration_test/conftest.py
+++ b/integration_test/conftest.py
@@ -26,6 +26,12 @@ def yatiml_log_warning():
     yatiml.logger.setLevel(logging.WARNING)
 
 
+def ls_snapshots(run_dir, instance=None):
+    """List all snapshots of the instance or workflow"""
+    return sorted(run_dir.snapshot_dir(instance).iterdir(),
+                  key=lambda path: tuple(map(int, path.stem.split("_")[1:])))
+
+
 def start_mmp_server(control_pipe, ymmsl_doc, run_dir):
     control_pipe[0].close()
     manager = Manager(ymmsl_doc, run_dir)
@@ -108,7 +114,8 @@ def run_manager_with_actors(
         for instance_name, callable in python_actors.items():
             proc = mp.Process(
                     target=_python_wrapper,
-                    args=(instance_name, env['MUSCLE_MANAGER'], callable))
+                    args=(instance_name, env['MUSCLE_MANAGER'], callable),
+                    name=instance_name)
             proc.start()
             python_processes.append(proc)
 
diff --git a/integration_test/test_cpp_mpp_client.py b/integration_test/test_cpp_mpp_client.py
index a08e7bc4..976dc41c 100644
--- a/integration_test/test_cpp_mpp_client.py
+++ b/integration_test/test_cpp_mpp_client.py
@@ -23,7 +23,7 @@ def tcp_server_process(control_pipe):
     message = MPPMessage(
             Reference('test_sender.test_port'),
             receiver,
-            10, 1.0, 2.0, settings, data).encoded()
+            10, 1.0, 2.0, settings, 0, 1.0, data).encoded()
 
     def handle_request(request_bytes):
         request = msgpack.unpackb(request_bytes, raw=False)
diff --git a/integration_test/test_snapshot_complex_coupling.py b/integration_test/test_snapshot_complex_coupling.py
new file mode 100644
index 00000000..dad2ee34
--- /dev/null
+++ b/integration_test/test_snapshot_complex_coupling.py
@@ -0,0 +1,188 @@
+import random
+import time
+
+import pytest
+from ymmsl import KeepsStateForNextUse, Operator, load, dump
+
+from libmuscle import Instance, Message
+from libmuscle.manager.run_dir import RunDir
+
+from .conftest import run_manager_with_actors, ls_snapshots
+
+
+_LOG_LEVEL = 'INFO'  # set to DEBUG for additional debug info
+
+
+def cache_component(max_channels=2):
+    ports = {Operator.F_INIT: [f'in{i+1}' for i in range(max_channels)],
+             Operator.O_I: [f'sub_out{i+1}' for i in range(max_channels)],
+             Operator.S: [f'sub_in{i+1}' for i in range(max_channels)],
+             Operator.O_F: [f'out{i+1}' for i in range(max_channels)]}
+    instance = Instance(ports)
+
+    cache_t = float('-inf')
+    cache_data = []
+    max_cache_age = None
+    nil_msg = Message(0.0, None, None)
+
+    while instance.reuse_instance():
+        if instance.resuming():
+            instance.load_snapshot()
+
+        if instance.should_init():
+            cache_valid_range = instance.get_setting('cache_valid', '[float]')
+            if max_cache_age is None:
+                max_cache_age = random.uniform(*cache_valid_range)
+
+            msgs = [instance.receive(port, default=nil_msg)
+                    for port in ports[Operator.F_INIT]]
+            cur_t = msgs[0].timestamp
+
+        if cur_t - cache_t >= max_cache_age:
+            # Cached value is no longer valid, run submodel for updated data
+            for msg, port in zip(msgs, ports[Operator.O_I]):
+                instance.send(port, Message(cur_t, None, msg.data))
+            cache_data = [instance.receive(port, default=nil_msg).data
+                          for port in ports[Operator.S]]
+            cache_t = cur_t
+            max_cache_age = random.uniform(*cache_valid_range)
+
+        for data, port in zip(cache_data, ports[Operator.O_F]):
+            instance.send(port, Message(cur_t, None, data))
+
+        if instance.should_save_final_snapshot():
+            instance.save_final_snapshot(Message(cur_t, None, []))
+
+
+def echo_component(max_channels=2):
+    ports = {Operator.F_INIT: [f'in{i+1}' for i in range(max_channels)],
+             Operator.O_F: [f'out{i+1}' for i in range(max_channels)]}
+    instance = Instance(ports, keeps_state_for_next_use=KeepsStateForNextUse.NO)
+
+    while instance.reuse_instance():
+        for p_in, p_out in zip(ports[Operator.F_INIT], ports[Operator.O_F]):
+            if instance.is_connected(p_in):
+                instance.send(p_out, instance.receive(p_in))
+
+
+def main_component():
+    instance = Instance({
+            Operator.O_I: ['state_out'],
+            Operator.S: ['Ai', 'Bi', 'Ci', 'Di'],
+            Operator.O_F: ['o_f']})
+
+    while instance.reuse_instance():
+        dt = instance.get_setting('dt', 'float')
+        t_max = instance.get_setting('t_max', 'float')
+
+        if instance.resuming():
+            msg = instance.load_snapshot()
+            t_cur = msg.timestamp
+            i, t_remaining = msg.data
+            monotonic_end = time.monotonic() + t_remaining
+
+        if instance.should_init():
+            t_cur = 0
+            monotonic_end = time.monotonic() + t_max
+            i = 0
+
+        while time.monotonic() < monotonic_end:
+            instance.send('state_out', Message(t_cur, None, i))
+            for port in ('Ai', 'Bi', 'Ci', 'Di'):
+                instance.receive(port)
+
+            t_cur += dt
+            i += 1
+            time.sleep(0.05)
+
+            if instance.should_save_snapshot(t_cur):
+                instance.save_snapshot(Message(
+                        t_cur, None, [i, monotonic_end - time.monotonic()]))
+
+        instance.send('o_f', Message(t_cur, None, i))
+
+        if instance.should_save_final_snapshot():
+            instance.save_final_snapshot(Message(t_cur, None, [i, 0]))
+
+
+@pytest.fixture
+def config():
+    return load(f"""ymmsl_version: v0.1
+model:
+  name: test_snapshot
+  components:
+    main: main_component
+    cacheA: cache_component
+    cacheB: cache_component
+    cacheC: cache_component
+    calcA: echo_component
+    calcB: echo_component
+    calcC: echo_component
+    calcD: echo_component
+  conduits:
+    main.state_out:
+    - cacheA.in1
+    - cacheB.in1
+    - cacheC.in1
+    - calcD.in1
+
+    cacheA.out1: main.Ai
+    cacheA.out2: main.Bi
+    cacheA.sub_out1: calcA.in1
+    cacheA.sub_out2: calcA.in2
+    calcA.out1: cacheA.sub_in1
+    calcA.out2: cacheA.sub_in2
+
+    cacheB.out1:
+    - cacheA.in2
+    - cacheC.in2
+    cacheB.sub_out1: calcB.in1
+    calcB.out1: cacheB.sub_in1
+
+    cacheC.out1: main.Ci
+    cacheC.sub_out1: calcC.in1
+    cacheC.sub_out2: calcC.in2
+    calcC.out1: cacheC.sub_in1
+
+    calcD.out1: main.Di
+
+settings:
+  dt: 1.234
+  t_max: 2.0  # seconds
+  cacheA.cache_valid: [2.0, 5.0]
+  cacheB.cache_valid: [3.0, 8.0]
+  cacheC.cache_valid: [4.0, 10.0]
+  muscle_remote_log_level: {_LOG_LEVEL}
+checkpoints:
+  at_end: true
+  wallclock_time:
+  - every: 0.5""")
+
+
+def test_snapshot_complex_coupling(tmp_path, config):
+    actors = {'main': main_component}
+    for c in 'ABC':
+        actors['cache' + c] = cache_component
+    for c in 'ABCD':
+        actors['calc' + c] = echo_component
+
+    run_dir1 = RunDir(tmp_path / 'run1')
+    run_manager_with_actors(
+            dump(config), run_dir1.path, python_actors=actors)
+
+    assert len(ls_snapshots(run_dir1, 'main')) == 5  # 2.0/0.5, at_end
+    assert len(ls_snapshots(run_dir1, 'cacheA')) == 5  # 2.0/0.5, at_end
+    assert len(ls_snapshots(run_dir1, 'cacheB')) == 5  # 2.0/0.5, at_end
+    assert len(ls_snapshots(run_dir1, 'cacheC')) == 5  # 2.0/0.5, at_end
+    # Due to caches, calcA/B/C may not run every 0.5 seconds
+    assert 1 <= len(ls_snapshots(run_dir1, 'calcA')) <= 5
+    assert 1 <= len(ls_snapshots(run_dir1, 'calcB')) <= 5
+    assert 1 <= len(ls_snapshots(run_dir1, 'calcC')) <= 5
+    assert len(ls_snapshots(run_dir1, 'calcD')) == 5  # 2.0/0.5, at_end
+
+    snapshots_ymmsl = ls_snapshots(run_dir1)
+    snapshot_docs = list(map(load, snapshots_ymmsl))
+    # Snapshots based on wallclock time are less reliable. There is at least one
+    # resume yMMSL: the at_end collection. At most 4 more, one for each
+    # wallclock_time checkpoint.
+    assert 1 <= len(snapshot_docs) <= 5
diff --git a/integration_test/test_snapshot_dispatch.py b/integration_test/test_snapshot_dispatch.py
new file mode 100644
index 00000000..7102a43c
--- /dev/null
+++ b/integration_test/test_snapshot_dispatch.py
@@ -0,0 +1,131 @@
+import pytest
+from ymmsl import KeepsStateForNextUse, Operator, load, dump
+
+from libmuscle import Instance, Message
+from libmuscle.manager.run_dir import RunDir
+
+from .conftest import run_manager_with_actors, ls_snapshots
+
+
+_LOG_LEVEL = 'INFO'  # set to DEBUG for additional debug info
+
+
+def component():
+    instance = Instance({
+            Operator.F_INIT: ['f_i'],
+            Operator.O_F: ['o_f']})
+
+    while instance.reuse_instance():
+        dt = instance.get_setting('dt', 'float')
+        t_max = instance.get_setting('t_max', 'float')
+
+        if instance.resuming():
+            msg = instance.load_snapshot()
+            t_cur = msg.timestamp
+            i, t_stop = msg.data
+
+        if instance.should_init():
+            msg = instance.receive('f_i', default=Message(0, None, 0))
+            t_cur = msg.timestamp
+            i = msg.data
+            t_stop = t_cur + t_max
+
+        while t_cur < t_stop:
+            # faux time-integration for testing snapshots
+            t_cur += dt
+
+            if instance.should_save_snapshot(t_cur):
+                instance.save_snapshot(Message(t_cur, None, [i, t_stop]))
+
+        instance.send('o_f', Message(t_cur, None, i))
+
+        if instance.should_save_final_snapshot():
+            instance.save_final_snapshot(Message(t_cur, None, [i, t_stop]))
+
+
+def stateless_component():
+    instance = Instance({
+            Operator.F_INIT: ['f_i'],
+            Operator.O_F: ['o_f']},
+            stateful=KeepsStateForNextUse.NO)
+
+    while instance.reuse_instance():
+        dt = instance.get_setting('dt', 'float')
+        t_max = instance.get_setting('t_max', 'float')
+
+        msg = instance.receive('f_i', default=Message(0, None, 0))
+        t_cur = msg.timestamp
+        i = msg.data
+        t_stop = t_cur + t_max
+
+        while t_cur < t_stop:
+            # faux time-integration for testing snapshots
+            t_cur += dt
+
+        instance.send('o_f', Message(t_cur, None, i))
+
+
+@pytest.fixture
+def dispatch_config():
+    return load(f"""ymmsl_version: v0.1
+model:
+  name: test_snapshot
+  components:
+    comp1: component
+    comp2: component
+    comp3: component
+    comp4: component
+    comp5: component
+  conduits:
+    comp1.o_f: comp2.f_i
+    comp2.o_f: comp3.f_i
+    comp3.o_f: comp4.f_i
+    comp4.o_f: comp5.f_i
+settings:
+  dt: 0.1234
+  t_max: 2.0
+  muscle_remote_log_level: {_LOG_LEVEL}
+checkpoints:
+  at_end: true
+  simulation_time:
+  - every: 2.5
+  - at:
+    - 2.3
+    - 2.8""")
+
+
+def test_snapshot_dispatch(tmp_path, dispatch_config):
+    actors = {f'comp{i + 1}': component for i in range(5)}
+    run_dir1 = RunDir(tmp_path / 'run1')
+    run_manager_with_actors(
+            dump(dispatch_config), run_dir1.path, python_actors=actors)
+
+    assert len(ls_snapshots(run_dir1, 'comp1')) == 2  # t=0, at_end
+    assert len(ls_snapshots(run_dir1, 'comp2')) == 5  # t=0, 2.5, 2.3, 2.8, at_end
+    assert len(ls_snapshots(run_dir1, 'comp3')) == 3  # t=2.5, 5, at_end
+    assert len(ls_snapshots(run_dir1, 'comp4')) == 3  # t=5, 7.5, at_end
+    assert len(ls_snapshots(run_dir1, 'comp5')) == 3  # t=7.5, 10, at_end
+
+    snapshots_ymmsl = ls_snapshots(run_dir1)
+    snapshot_docs = list(map(load, snapshots_ymmsl))
+    # More ymmsl restarts files may be possible, depending on the sequence of
+    # incoming SnapshotMetadata...
+    assert len(snapshot_docs) == 16
+
+    # resume from the snapshots taken at t>=2.3
+    run_dir2 = RunDir(tmp_path / 'run2')
+    dispatch_config.resume = {
+        'comp1': ls_snapshots(run_dir1, 'comp1')[1],
+        'comp2': ls_snapshots(run_dir1, 'comp2')[1]}
+
+    run_manager_with_actors(
+            dump(dispatch_config), run_dir2.path, python_actors=actors)
+
+    assert len(ls_snapshots(run_dir2, 'comp1')) == 1  # resume
+    assert len(ls_snapshots(run_dir2, 'comp2')) == 4  # resume, t=2.5, 2.8, at_end
+    assert len(ls_snapshots(run_dir2, 'comp3')) == 3  # t=2.5, 5, at_end
+    assert len(ls_snapshots(run_dir2, 'comp4')) == 3  # t=5, 7.5, at_end
+    assert len(ls_snapshots(run_dir2, 'comp5')) == 3  # t=7.5, 10, at_end
+    # More ymmsl restarts files may be possible, depending on the sequence of
+    # incoming SnapshotMetadata...
+    assert len(ls_snapshots(run_dir2)) >= 13
diff --git a/integration_test/test_snapshot_interact.py b/integration_test/test_snapshot_interact.py
new file mode 100644
index 00000000..5492f9e2
--- /dev/null
+++ b/integration_test/test_snapshot_interact.py
@@ -0,0 +1,164 @@
+import logging
+import sys
+from pathlib import Path
+
+import pytest
+from ymmsl import Operator, load, dump
+
+from libmuscle import Instance, Message
+from libmuscle.manager.run_dir import RunDir
+
+from .conftest import run_manager_with_actors, ls_snapshots
+
+# Make interact_coupling.py available (from docs/sources/examples)
+sys.path.append(str(
+        Path(__file__).parents[1] / 'docs' / 'source' / 'examples' / 'python'))
+import interact_coupling  # noqa
+
+_LOG_LEVEL = 'INFO'  # set to DEBUG for additional debug info
+
+
+def component():
+    instance = Instance({
+            Operator.O_I: ['o_i'],
+            Operator.S: ['s']})
+
+    while instance.reuse_instance():
+        t0 = instance.get_setting('t0', 'float')
+        dt = instance.get_setting('dt', 'float')
+        t_max = instance.get_setting('t_max', 'float')
+
+        if instance.resuming():
+            msg = instance.load_snapshot()
+            t_cur = msg.timestamp
+            i, t_stop = msg.data
+
+        if instance.should_init():
+            t_cur = t0
+            i = 0
+            t_stop = t0 + t_max
+
+        rcvd_i = 0
+        while t_cur < t_stop:
+            # faux time-integration for testing snapshots
+            t_next = t_cur + dt
+            if t_next >= t_stop:
+                t_next = None
+            logging.info(f'Sending {i} at {t_cur}, next at {t_next}')
+            instance.send('o_i', Message(t_cur, t_next, i))
+
+            msg = instance.receive('s')
+            logging.info(
+                    f'Received {msg.data} from time {msg.timestamp},'
+                    f' next at {msg.next_timestamp}')
+            assert msg.data >= rcvd_i
+            rcvd_i = msg.data
+
+            t_cur += dt
+            i += 1
+
+            if instance.should_save_snapshot(t_cur):
+                instance.save_snapshot(Message(t_cur, None, [i, t_stop]))
+
+        if instance.should_save_final_snapshot():
+            instance.save_final_snapshot(Message(t_cur, None, [i, t_stop]))
+
+
+def test_snapshot_interact_lockstep(tmp_path):
+    config = f"""ymmsl_version: v0.1
+model:
+  name: test_snapshot
+  components:
+    comp1: component
+    comp2: component
+  conduits:
+    comp1.o_i: comp2.s
+    comp2.o_i: comp1.s
+settings:
+  t0: 0.35
+  dt: 0.1234
+  t_max: 3.0
+  muscle_remote_log_level: {_LOG_LEVEL}
+checkpoints:
+  simulation_time:
+  - every: 1.0
+    start: 0.75
+    stop: 2.0
+  - at:
+    - 2.5"""
+    actors = {f'comp{i + 1}': component for i in range(2)}
+
+    run_dir1 = RunDir(tmp_path / 'run1')
+    run_manager_with_actors(config, run_dir1.path, python_actors=actors)
+
+    assert len(ls_snapshots(run_dir1, 'comp1')) == 3  # t=0.75, 1.75, 2.5
+    assert len(ls_snapshots(run_dir1, 'comp2')) == 3  # t=0.75, 1.75, 2.5
+
+    snapshots_ymmsl = ls_snapshots(run_dir1)
+    snapshot_docs = list(map(load, snapshots_ymmsl))
+    assert len(snapshot_docs) == 3
+
+    # resume from the snapshots taken at t>=1.75
+    run_dir2 = RunDir(tmp_path / 'run2')
+    config_doc = load(config)
+    config_doc.update(snapshot_docs[1])
+
+    run_manager_with_actors(
+            dump(config_doc), run_dir2.path, python_actors=actors)
+
+    assert len(ls_snapshots(run_dir2, 'comp1')) == 2  # resume, t=2.5
+    assert len(ls_snapshots(run_dir2, 'comp2')) == 2  # resume, t=2.5
+    assert len(ls_snapshots(run_dir2)) == 2
+
+
+@pytest.mark.parametrize('scale', [0.1, 0.9, 1.0, 1.1, 1.5])
+def test_snapshot_interact_varstep(tmp_path, scale):
+    config = f"""ymmsl_version: v0.1
+model:
+  name: test_snapshot
+  components:
+    comp1: component
+    comp2: component
+    coupler: checkpointing_temporal_coupler
+  conduits:
+    comp1.o_i: coupler.a_in
+    coupler.a_out: comp1.s
+    comp2.o_i: coupler.b_in
+    coupler.b_out: comp2.s
+settings:
+  t0: 0.35
+  comp1.dt: 0.1234
+  comp2.dt: {0.1234 * scale}
+  t_max: 3.0
+  muscle_remote_log_level: {_LOG_LEVEL}
+checkpoints:
+  simulation_time:
+  - every: 1.0
+    start: 0.75
+    stop: 2.0
+  - at:
+    - 2.5"""
+    actors = {f'comp{i + 1}': component for i in range(2)}
+    actors['coupler'] = interact_coupling.checkpointing_temporal_coupler
+
+    run_dir1 = RunDir(tmp_path / 'run1')
+    run_manager_with_actors(config, run_dir1.path, python_actors=actors)
+
+    assert len(ls_snapshots(run_dir1, 'comp1')) == 3  # t=0.75, 1.75, 2.5
+    assert len(ls_snapshots(run_dir1, 'comp2')) == 3  # t=0.75, 1.75, 2.5
+
+    snapshots_ymmsl = ls_snapshots(run_dir1)
+    snapshot_docs = list(map(load, snapshots_ymmsl))
+    assert len(snapshot_docs) == 3
+
+    # resume from the snapshots taken at t>=1.75
+    run_dir2 = RunDir(tmp_path / 'run2')
+    config_doc = load(config)
+    config_doc.update(snapshot_docs[1])
+
+    run_manager_with_actors(
+            dump(config_doc), run_dir2.path, python_actors=actors)
+
+    assert len(ls_snapshots(run_dir2, 'comp1')) == 2  # resume, t=2.5
+    assert len(ls_snapshots(run_dir2, 'comp2')) == 2  # resume, t=2.5
+    assert len(ls_snapshots(run_dir2)) == 2
diff --git a/integration_test/test_snapshot_macro_micro.py b/integration_test/test_snapshot_macro_micro.py
new file mode 100644
index 00000000..885ac704
--- /dev/null
+++ b/integration_test/test_snapshot_macro_micro.py
@@ -0,0 +1,305 @@
+import pytest
+from ymmsl import KeepsStateForNextUse, Operator, load, dump
+
+from libmuscle import Instance, Message
+from libmuscle.manager.run_dir import RunDir
+
+from .conftest import run_manager_with_actors, ls_snapshots
+
+
+_LOG_LEVEL = 'INFO'  # set to DEBUG for additional debug info
+
+
+def macro():
+    instance = Instance({
+            Operator.O_I: ['o_i'],
+            Operator.S: ['s']})
+
+    while instance.reuse_instance():
+        dt = instance.get_setting('dt', 'float')
+        t_max = instance.get_setting('t_max', 'float')
+
+        if instance.resuming():
+            msg = instance.load_snapshot()
+            # load state from message
+            t_cur = msg.timestamp
+            i = msg.data
+            assert i >= 1
+
+        if instance.should_init():
+            t_cur = instance.get_setting('t0', 'float')
+            i = 0
+
+        while t_cur + dt <= t_max:
+            t_next = t_cur + dt
+            if t_next + dt > t_max:
+                t_next = None  # final iteration of this time-integration loop
+            instance.send('o_i', Message(t_cur, t_next, i))
+
+            msg = instance.receive('s')
+            assert msg.data == i
+
+            i += 1
+            t_cur += dt
+
+            if instance.should_save_snapshot(t_cur):
+                instance.save_snapshot(Message(t_cur, None, i))
+
+        if instance.should_save_final_snapshot():
+            instance.save_final_snapshot(Message(t_cur, None, i))
+
+
+def macro_vector():
+    instance = Instance({
+            Operator.O_I: ['o_i[]'],
+            Operator.S: ['s[]']})
+
+    while instance.reuse_instance():
+        dt = instance.get_setting('dt', 'float')
+        t_max = instance.get_setting('t_max', 'float')
+
+        if instance.resuming():
+            msg = instance.load_snapshot()
+            # load state from message
+            t_cur = msg.timestamp
+            i = msg.data
+            assert i >= 1
+
+        if instance.should_init():
+            t_cur = instance.get_setting('t0', 'float')
+            i = 0
+
+        while t_cur + dt <= t_max:
+            t_next = t_cur + dt
+            if t_next + dt > t_max:
+                t_next = None  # final iteration of this time-integration loop
+            for slot in range(instance.get_port_length('o_i')):
+                instance.send('o_i', Message(t_cur, t_next, i), slot)
+
+            for slot in range(instance.get_port_length('s')):
+                msg = instance.receive('s', slot)
+                assert msg.data == i
+
+            i += 1
+            t_cur += dt
+
+            if instance.should_save_snapshot(t_cur):
+                instance.save_snapshot(Message(t_cur, None, i))
+
+        if instance.should_save_final_snapshot():
+            instance.save_final_snapshot(Message(t_cur, None, i))
+
+
+def micro():
+    instance = Instance({
+            Operator.F_INIT: ['f_i'],
+            Operator.O_F: ['o_f']})
+
+    while instance.reuse_instance():
+        dt = instance.get_setting('dt', 'float')
+        t_max = instance.get_setting('t_max', 'float')
+
+        if instance.resuming():
+            msg = instance.load_snapshot()
+            t_cur = msg.timestamp
+            i, t_stop = msg.data
+
+        if instance.should_init():
+            msg = instance.receive('f_i')
+            t_cur = msg.timestamp
+            i = msg.data
+            t_stop = t_cur + t_max
+
+        while t_cur < t_stop:
+            # faux time-integration for testing snapshots
+            t_cur += dt
+
+            if instance.should_save_snapshot(t_cur):
+                instance.save_snapshot(Message(t_cur, None, [i, t_stop]))
+
+        instance.send('o_f', Message(t_cur, None, i))
+
+        if instance.should_save_final_snapshot():
+            instance.save_final_snapshot(Message(t_cur, None, [i, t_stop]))
+
+
+def stateless_micro():
+    instance = Instance({
+            Operator.F_INIT: ['f_i'],
+            Operator.O_F: ['o_f']},
+            keeps_state_for_next_use=KeepsStateForNextUse.NO)
+
+    while instance.reuse_instance():
+        dt = instance.get_setting('dt', 'float')
+        t_max = instance.get_setting('t_max', 'float')
+
+        msg = instance.receive('f_i')
+        t_cur = msg.timestamp
+        i = msg.data
+        t_stop = t_cur + t_max
+
+        while t_cur < t_stop:
+            # faux time-integration for testing snapshots
+            t_cur += dt
+
+        instance.send('o_f', Message(t_cur, None, i))
+
+
+def data_transformer():
+    instance = Instance({
+            Operator.F_INIT: ['f_i'],
+            Operator.O_F: ['o_f']},
+            keeps_state_for_next_use=KeepsStateForNextUse.NO)
+
+    while instance.reuse_instance():
+        msg = instance.receive('f_i')
+        instance.send('o_f', msg)
+
+
+@pytest.fixture
+def base_config():
+    return load(f"""ymmsl_version: v0.1
+model:
+  name: test_snapshot
+  components:
+    macro: macro_implementation
+    micro: micro_implementation
+  conduits:
+    macro.o_i: micro.f_i
+    micro.o_f: macro.s
+settings:
+  macro.t0: 0.14
+  macro.dt: 0.17
+  macro.t_max: 1.9
+  micro.dt: 0.009
+  micro.t_max: 0.1
+  muscle_remote_log_level: {_LOG_LEVEL}
+checkpoints:
+  at_end: true
+  simulation_time:
+  - every: 0.4""")
+
+
+@pytest.fixture
+def config_with_transformer(base_config):
+    base_config.update(load("""ymmsl_version: v0.1
+model:
+  name: test_snapshot
+  components:
+    transformer1: transformer
+    transformer2: transformer
+  conduits:
+    macro.o_i: transformer1.f_i
+    transformer1.o_f: micro.f_i
+    micro.o_f: transformer2.f_i
+    transformer2.o_f: macro.s"""))
+    return base_config
+
+
+def test_snapshot_macro_micro(tmp_path, base_config):
+    actors = {'macro': macro, 'micro': micro}
+    run_dir1 = RunDir(tmp_path / 'run1')
+    run_manager_with_actors(
+            dump(base_config), run_dir1.path, python_actors=actors)
+
+    macro_snapshots = ls_snapshots(run_dir1, 'macro')
+    assert len(macro_snapshots) == 6  # 0, 0.4, 0.8, 1.2, 1.6, final
+    micro_snapshots = ls_snapshots(run_dir1, 'micro')
+    assert len(micro_snapshots) == 6  # 0, 0.4, 0.8, 1.2, 1.6, final
+    snapshots_ymmsl = ls_snapshots(run_dir1)
+    snapshot_docs = list(map(load, snapshots_ymmsl))
+    assert len(snapshot_docs) == 7
+    assert snapshot_docs[0].resume['macro'] == macro_snapshots[0]
+    assert snapshot_docs[0].resume['micro'] == micro_snapshots[0]
+    assert snapshot_docs[1].resume['macro'] == macro_snapshots[0]
+    assert snapshot_docs[1].resume['micro'] == micro_snapshots[1]
+    for i in range(2, 7):
+        assert snapshot_docs[i].resume['macro'] == macro_snapshots[i - 1]
+        assert snapshot_docs[i].resume['micro'] == micro_snapshots[i - 1]
+
+    # resume from the snapshots taken at t>=1.2
+    run_dir2 = RunDir(tmp_path / 'run2')
+    base_config.update(snapshot_docs[4])  # add resume info
+    run_manager_with_actors(
+            dump(base_config), run_dir2.path, python_actors=actors)
+
+    assert len(ls_snapshots(run_dir2, 'macro')) == 3  # resume, 1.6, final
+    assert len(ls_snapshots(run_dir2, 'micro')) == 3  # resume, 1.6, final
+    assert len(ls_snapshots(run_dir2)) == 3
+
+    # resume from the first workflow snapshot (this restarts macro from scratch)
+    run_dir3 = RunDir(tmp_path / 'run3')
+    base_config.resume = {}                     # clear resume information
+    base_config.update(snapshot_docs[0])        # add resume info
+    base_config.settings['macro.t_max'] = 0.6   # run shorter
+    run_manager_with_actors(
+            dump(base_config), run_dir3.path, python_actors=actors)
+
+
+def test_snapshot_macro_stateless_micro(tmp_path, base_config):
+    actors = {'macro': macro, 'micro': stateless_micro}
+    run_dir1 = RunDir(tmp_path / 'run1')
+    run_manager_with_actors(
+            dump(base_config), run_dir1.path, python_actors=actors)
+
+    assert len(ls_snapshots(run_dir1, 'macro')) == 6  # 0, 0.4, 0.8, 1.2, 1.6, final
+    assert len(ls_snapshots(run_dir1, 'micro')) == 6  # 0, 0.4, 0.8, 1.2, 1.6, final
+    snapshots_ymmsl = ls_snapshots(run_dir1)
+    snapshot_docs = list(map(load, snapshots_ymmsl))
+    assert len(snapshot_docs) == 6
+
+    # resume from the snapshot taken at t>=1.2
+    run_dir2 = RunDir(tmp_path / 'run2')
+    base_config.update(snapshot_docs[3])  # add resume info
+    run_manager_with_actors(
+            dump(base_config), run_dir2.path, python_actors=actors)
+
+    assert len(ls_snapshots(run_dir2, 'macro')) == 3  # resume, 1.6, final
+    assert len(ls_snapshots(run_dir2, 'micro')) == 4  # resume, 1.2, 1.6, final
+    assert len(ls_snapshots(run_dir2)) == 3
+
+
+def test_snapshot_macro_vector_micro(tmp_path, base_config):
+    base_config.model.components[1].multiplicity = [2]
+    actors = {'macro': macro_vector, 'micro[0]': micro, 'micro[1]': micro}
+
+    run_dir1 = RunDir(tmp_path / 'run1')
+    run_manager_with_actors(
+            dump(base_config), run_dir1.path, python_actors=actors)
+
+    assert len(ls_snapshots(run_dir1, 'macro')) == 6  # 0, 0.4, 0.8, 1.2, 1.6, final
+    assert len(ls_snapshots(run_dir1, 'micro[0]')) == 6  # 0, 0.4, 0.8, 1.2, 1.6, final
+    assert len(ls_snapshots(run_dir1, 'micro[1]')) == 6  # 0, 0.4, 0.8, 1.2, 1.6, final
+    snapshots_ymmsl = ls_snapshots(run_dir1)
+    assert len(snapshots_ymmsl) == 8
+
+    run_dir2 = RunDir(tmp_path / 'run2')
+    base_config.update(load(snapshots_ymmsl[-3]))  # add resume info
+    run_manager_with_actors(
+            dump(base_config), run_dir2.path, python_actors=actors)
+
+    assert len(ls_snapshots(run_dir2, 'macro')) == 3  # resume, 1.6, final
+    assert len(ls_snapshots(run_dir2, 'micro[0]')) == 3  # resume, 1.6, final
+    assert len(ls_snapshots(run_dir2, 'micro[1]')) == 3  # resume, 1.6, final
+    assert len(ls_snapshots(run_dir2)) == 3
+
+
+def test_snapshot_macro_transformer_micro(tmp_path, config_with_transformer):
+    actors = {'macro': macro, 'micro': micro, 'transformer1': data_transformer,
+              'transformer2': data_transformer}
+
+    run_dir1 = RunDir(tmp_path / 'run1')
+    run_manager_with_actors(
+            dump(config_with_transformer), run_dir1.path, python_actors=actors)
+
+    snapshots_ymmsl = ls_snapshots(run_dir1)
+    assert len(snapshots_ymmsl) == 8
+
+    # pick one to resume from
+    run_dir2 = RunDir(tmp_path / 'run2')
+    config_with_transformer.update(load(snapshots_ymmsl[4]))  # add resume info
+    run_manager_with_actors(
+            dump(config_with_transformer), run_dir2.path, python_actors=actors)
+
+    snapshots_ymmsl = ls_snapshots(run_dir2)
+    assert len(snapshots_ymmsl) == 6
diff --git a/libmuscle/cpp/src/libmuscle/tests/mocks/mock_mpp_client.cpp b/libmuscle/cpp/src/libmuscle/tests/mocks/mock_mpp_client.cpp
index c23dbafc..31507a44 100644
--- a/libmuscle/cpp/src/libmuscle/tests/mocks/mock_mpp_client.cpp
+++ b/libmuscle/cpp/src/libmuscle/tests/mocks/mock_mpp_client.cpp
@@ -44,7 +44,8 @@ Settings MockMPPClient::make_overlay_() {
 }
 
 MPPMessage MockMPPClient::next_receive_message(
-        "test.out", "test2.in", 0, 0.0, 1.0, make_overlay_(), Data::dict("test1", 12));
+        "test.out", "test2.in", 0, 0.0, 1.0, make_overlay_(),0, 9.0,
+        Data::dict("test1", 12));
 
 Reference MockMPPClient::last_receiver("_none");
 
diff --git a/libmuscle/cpp/src/libmuscle/tests/mocks/mock_post_office.cpp b/libmuscle/cpp/src/libmuscle/tests/mocks/mock_post_office.cpp
index e1d66eac..cb6a92d5 100644
--- a/libmuscle/cpp/src/libmuscle/tests/mocks/mock_post_office.cpp
+++ b/libmuscle/cpp/src/libmuscle/tests/mocks/mock_post_office.cpp
@@ -6,13 +6,13 @@ int MockPostOffice::handle_request(
         char const * res_buf, std::size_t res_len,
         std::unique_ptr<DataConstRef> & response) {
     response = std::make_unique<DataConstRef>(
-            MPPMessage("test.out", "test2.in", 0, 0.0, 1.0, Data(), Data()).encoded());
+            MPPMessage("test.out", "test2.in", 0, 0.0, 1.0, Data(), 0, 0.0, Data()).encoded());
     return -1;
 }
 
 std::unique_ptr<DataConstRef> MockPostOffice::get_response(int fd) {
     return std::make_unique<DataConstRef>(
-            MPPMessage("test.out", "test2.in", 0, 0.0, 1.0, Data(), Data()).encoded());
+            MPPMessage("test.out", "test2.in", 0, 0.0, 1.0, Data(), 0, 8.0, Data()).encoded());
 }
 
 void MockPostOffice::deposit(
diff --git a/libmuscle/cpp/src/libmuscle/tests/tcp_transport_server_test.cpp b/libmuscle/cpp/src/libmuscle/tests/tcp_transport_server_test.cpp
index 9e069031..7082a0bc 100644
--- a/libmuscle/cpp/src/libmuscle/tests/tcp_transport_server_test.cpp
+++ b/libmuscle/cpp/src/libmuscle/tests/tcp_transport_server_test.cpp
@@ -38,6 +38,7 @@ int main(int argc, char *argv[]) {
             "test_sender.port", receiver, 10,
             0.0, 1.0,
             overlay_settings,
+            0, 6.0,
             data_dict);
     auto msg_data = std::make_unique<DataConstRef>(msg.encoded());
     post_office.deposit(receiver, std::move(msg_data));
diff --git a/libmuscle/cpp/src/libmuscle/tests/test_mpp_message.cpp b/libmuscle/cpp/src/libmuscle/tests/test_mpp_message.cpp
index 88ba96c9..cebcedd5 100644
--- a/libmuscle/cpp/src/libmuscle/tests/test_mpp_message.cpp
+++ b/libmuscle/cpp/src/libmuscle/tests/test_mpp_message.cpp
@@ -29,7 +29,7 @@ TEST(test_mcp_message, create_mcp_message) {
             Reference("sender.port"), Reference("receiver.port"),
             10,
             100.1, 101.0,
-            test, abc
+            test, 0, 1.0, abc
             );
 
     ASSERT_EQ(m.sender, "sender.port");
@@ -38,6 +38,8 @@ TEST(test_mcp_message, create_mcp_message) {
     ASSERT_EQ(m.timestamp, 100.1);
     ASSERT_EQ(m.next_timestamp, 101.0);
     ASSERT_EQ(m.settings_overlay.as<std::string>(), "test");
+    ASSERT_EQ(m.message_number, 0);
+    ASSERT_EQ(m.saved_until, 1.0);
     ASSERT_EQ(m.data.as<std::string>(), "abc");
 }
 
@@ -48,7 +50,7 @@ TEST(test_mcp_message, create_mcp_message_minimal) {
             Reference("sender.port"), Reference("receiver.port"),
             {},
             100.1, {},
-            test, abc
+            test, 0, 2.0, abc
             );
 
     ASSERT_EQ(m.sender, "sender.port");
@@ -57,6 +59,8 @@ TEST(test_mcp_message, create_mcp_message_minimal) {
     ASSERT_EQ(m.timestamp, 100.1);
     ASSERT_FALSE(m.next_timestamp.is_set());
     ASSERT_TRUE(m.settings_overlay.is_nil());
+    ASSERT_EQ(m.message_number, 0);
+    ASSERT_EQ(m.saved_until, 2.0);
     ASSERT_TRUE(m.data.is_nil());
 }
 
@@ -68,6 +72,8 @@ TEST(test_mcp_message, from_bytes) {
             "timestamp", 100.1,
             "next_timestamp", Data(),
             "settings_overlay", Data(),
+            "message_number", 0,
+            "saved_until", 3.0,
             "data", Data()
             );
 
@@ -84,6 +90,8 @@ TEST(test_mcp_message, from_bytes) {
     ASSERT_EQ(m.timestamp, 100.1);
     ASSERT_FALSE(m.next_timestamp.is_set());
     ASSERT_TRUE(m.settings_overlay.is_nil());
+    ASSERT_EQ(m.message_number, 0);
+    ASSERT_EQ(m.saved_until, 3.0);
     ASSERT_TRUE(m.data.is_nil());
 }
 
diff --git a/libmuscle/cpp/src/libmuscle/tests/test_outbox.cpp b/libmuscle/cpp/src/libmuscle/tests/test_outbox.cpp
index 07486cb3..e98bd423 100644
--- a/libmuscle/cpp/src/libmuscle/tests/test_outbox.cpp
+++ b/libmuscle/cpp/src/libmuscle/tests/test_outbox.cpp
@@ -36,6 +36,7 @@ TEST(libmuscle_outbox, test_deposit_retrieve_message) {
             Optional<int>(),
             0.0, 1.0,
             DataConstRef(),
+            0, 1.0,
             DataConstRef("testing"));
 
     auto message_data = std::make_unique<DataConstRef>(message.encoded());
diff --git a/libmuscle/cpp/src/libmuscle/tests/test_post_office.cpp b/libmuscle/cpp/src/libmuscle/tests/test_post_office.cpp
index 68af8bf7..bf6981d8 100644
--- a/libmuscle/cpp/src/libmuscle/tests/test_post_office.cpp
+++ b/libmuscle/cpp/src/libmuscle/tests/test_post_office.cpp
@@ -36,7 +36,7 @@ std::unique_ptr<DataConstRef> make_message() {
             "test_sender.port", "test_receiver.port",
             Optional<int>(),
             0.0, 1.0,
-            DataConstRef(), DataConstRef());
+            DataConstRef(), 0, 5.0, DataConstRef());
     return std::make_unique<DataConstRef>(msg.encoded());
 }
 
diff --git a/libmuscle/cpp/src/libmuscle/tests/test_tcp_communication.cpp b/libmuscle/cpp/src/libmuscle/tests/test_tcp_communication.cpp
index f9c60c30..2d152161 100644
--- a/libmuscle/cpp/src/libmuscle/tests/test_tcp_communication.cpp
+++ b/libmuscle/cpp/src/libmuscle/tests/test_tcp_communication.cpp
@@ -36,7 +36,7 @@ TEST(test_tcp_communication, send_receive) {
     MPPMessage msg(
             "test_sender.port", receiver, 10,
             0.0, 1.0,
-            Data::dict("par1", 13),
+            Data::dict("par1", 13), 1, 4.0,
             Data::dict("var1", 1, "var2", 2.0, "var3", "3"));
     auto msg_data = std::make_unique<DataConstRef>(msg.encoded());
     post_office.deposit(receiver, std::move(msg_data));
@@ -53,6 +53,8 @@ TEST(test_tcp_communication, send_receive) {
     ASSERT_EQ(m.timestamp, 0.0);
     ASSERT_EQ(m.next_timestamp, 1.0);
     ASSERT_EQ(m.settings_overlay["par1"].as<int>(), 13);
+    ASSERT_EQ(m.message_number, 1);
+    ASSERT_EQ(m.saved_until, 4.0);
     ASSERT_EQ(m.data["var1"].as<int>(), 1);
     ASSERT_EQ(m.data["var2"].as<double>(), 2.0);
     ASSERT_EQ(m.data["var3"].as<std::string>(), "3");
diff --git a/libmuscle/python/libmuscle/manager/test/conftest.py b/libmuscle/python/libmuscle/manager/test/conftest.py
index 9ba095dd..992a3950 100644
--- a/libmuscle/python/libmuscle/manager/test/conftest.py
+++ b/libmuscle/python/libmuscle/manager/test/conftest.py
@@ -1,12 +1,13 @@
 from pathlib import Path
 
 import pytest
-from ymmsl import (Component, Conduit, Configuration, Model, Reference,
-                   Settings)
+from ymmsl import Component, Conduit, Configuration, Model, Reference
 
 from libmuscle.manager.instance_registry import InstanceRegistry
 from libmuscle.manager.logger import Logger
 from libmuscle.manager.mmp_server import MMPRequestHandler
+from libmuscle.manager.run_dir import RunDir
+from libmuscle.manager.snapshot_registry import SnapshotRegistry
 from libmuscle.manager.topology_store import TopologyStore
 
 
@@ -18,18 +19,8 @@ def logger(tmpdir):
 
 
 @pytest.fixture
-def settings():
-    return Settings()
-
-
-@pytest.fixture
-def instance_registry():
-    return InstanceRegistry()
-
-
-@pytest.fixture
-def topology_store() -> TopologyStore:
-    config = Configuration(
+def mmp_configuration():
+    return Configuration(
             Model(
                 'test_model',
                 [
@@ -41,13 +32,29 @@ def topology_store() -> TopologyStore:
                     Conduit('micro.out', 'macro.in')
                 ]))
 
-    return TopologyStore(config)
+
+@pytest.fixture
+def instance_registry():
+    return InstanceRegistry()
 
 
 @pytest.fixture
-def mmp_request_handler(logger, settings, instance_registry, topology_store):
+def topology_store(mmp_configuration) -> TopologyStore:
+    return TopologyStore(mmp_configuration)
+
+
+@pytest.fixture
+def snapshot_registry(mmp_configuration, topology_store) -> SnapshotRegistry:
+    return SnapshotRegistry(mmp_configuration, None, topology_store)
+
+
+@pytest.fixture
+def mmp_request_handler(
+        logger, mmp_configuration, instance_registry, topology_store,
+        snapshot_registry):
     return MMPRequestHandler(
-            logger, settings, instance_registry, topology_store)
+            logger, mmp_configuration, instance_registry, topology_store,
+            snapshot_registry, None)
 
 
 @pytest.fixture
@@ -63,14 +70,16 @@ def loaded_instance_registry(instance_registry):
 
 @pytest.fixture
 def registered_mmp_request_handler(
-        logger, settings, loaded_instance_registry, topology_store):
+        logger, mmp_configuration, loaded_instance_registry, topology_store,
+        snapshot_registry):
     return MMPRequestHandler(
-            logger, settings, loaded_instance_registry, topology_store)
+            logger, mmp_configuration, loaded_instance_registry, topology_store,
+            snapshot_registry, None)
 
 
 @pytest.fixture
-def topology_store2() -> TopologyStore:
-    config = Configuration(
+def mmp_configuration2():
+    return Configuration(
             Model(
                 'test_model',
                 [
@@ -85,7 +94,15 @@ def topology_store2() -> TopologyStore:
                     Conduit('meso.out', 'macro.in')
                 ]))
 
-    return TopologyStore(config)
+
+@pytest.fixture
+def topology_store2(mmp_configuration2) -> TopologyStore:
+    return TopologyStore(mmp_configuration2)
+
+
+@pytest.fixture
+def snapshot_registry2(mmp_configuration2, topology_store) -> SnapshotRegistry:
+    return SnapshotRegistry(mmp_configuration2, None, topology_store)
 
 
 @pytest.fixture
@@ -109,6 +126,9 @@ def loaded_instance_registry2():
 
 @pytest.fixture
 def registered_mmp_request_handler2(
-        logger, settings, loaded_instance_registry2, topology_store2):
+        logger, mmp_configuration, loaded_instance_registry2, topology_store2,
+        snapshot_registry2, tmp_path):
     return MMPRequestHandler(
-            logger, settings, loaded_instance_registry2, topology_store2)
+            logger, mmp_configuration,
+            loaded_instance_registry2, topology_store2, snapshot_registry2,
+            RunDir(tmp_path))
diff --git a/libmuscle/python/libmuscle/manager/test/test_mmp_request_handler.py b/libmuscle/python/libmuscle/manager/test/test_mmp_request_handler.py
index 6c40e02e..bc61f0a0 100644
--- a/libmuscle/python/libmuscle/manager/test/test_mmp_request_handler.py
+++ b/libmuscle/python/libmuscle/manager/test/test_mmp_request_handler.py
@@ -1,14 +1,22 @@
+import dataclasses
+from pathlib import Path
+from unittest.mock import MagicMock
+
 import msgpack
-from ymmsl import Operator, Reference
+from ymmsl import (
+        Operator, Reference, Checkpoints, CheckpointRangeRule, CheckpointAtRule)
 
 from libmuscle.logging import LogLevel
 from libmuscle.manager.mmp_server import MMPRequestHandler
 from libmuscle.mcp.protocol import RequestType, ResponseType
+from libmuscle.snapshot import SnapshotMetadata
 
 
-def test_create_servicer(logger, settings, instance_registry,
-                         topology_store):
-    MMPRequestHandler(logger, settings, instance_registry, topology_store)
+def test_create_servicer(logger, mmp_configuration, instance_registry,
+                         topology_store, snapshot_registry):
+    MMPRequestHandler(
+            logger, mmp_configuration, instance_registry, topology_store,
+            snapshot_registry, None)
 
 
 def test_log_message(mmp_request_handler, caplog):
@@ -31,7 +39,7 @@ def test_log_message(mmp_request_handler, caplog):
     assert caplog.records[0].message == 'Testing log message'
 
 
-def test_get_settings(settings, mmp_request_handler):
+def test_get_settings(mmp_configuration, mmp_request_handler):
     request = [RequestType.GET_SETTINGS.value]
     encoded_request = msgpack.packb(request, use_bin_type=True)
 
@@ -42,12 +50,12 @@ def test_get_settings(settings, mmp_request_handler):
     assert decoded_result[0] == ResponseType.SUCCESS.value
     assert decoded_result[1] == {}
 
-    settings['test1'] = 13
-    settings['test2'] = 12.3
-    settings['test3'] = 'testing'
-    settings['test4'] = True
-    settings['test5'] = [2.3, 7.4]
-    settings['test6'] = [[1.0, 2.0], [2.0, 1.0]]
+    mmp_configuration.settings['test1'] = 13
+    mmp_configuration.settings['test2'] = 12.3
+    mmp_configuration.settings['test3'] = 'testing'
+    mmp_configuration.settings['test4'] = True
+    mmp_configuration.settings['test5'] = [2.3, 7.4]
+    mmp_configuration.settings['test6'] = [[1.0, 2.0], [2.0, 1.0]]
 
     result = mmp_request_handler.handle_request(encoded_request)
     decoded_result = msgpack.unpackb(result, raw=False)
@@ -63,7 +71,7 @@ def test_get_settings(settings, mmp_request_handler):
     assert result_dict['test4'] is True
     assert result_dict['test5'] == [2.3, 7.4]
     assert result_dict['test6'] == [[1.0, 2.0], [2.0, 1.0]]
-    assert result_dict == settings.as_ordered_dict()
+    assert result_dict == mmp_configuration.settings.as_ordered_dict()
 
 
 def test_register_instance(mmp_request_handler, instance_registry):
@@ -86,6 +94,52 @@ def test_register_instance(mmp_request_handler, instance_registry):
     assert registered_ports['test_instance'][0].operator == Operator.F_INIT
 
 
+def test_get_checkpoint_info(mmp_configuration, mmp_request_handler):
+    resume_path = Path('/path/to/resume.pack')
+    mmp_configuration.resume = {Reference('test_instance'): resume_path}
+    mmp_configuration.checkpoints = Checkpoints(
+            True,
+            [CheckpointRangeRule(every=10), CheckpointAtRule([1, 2, 3.0])])
+
+    request = [RequestType.GET_CHECKPOINT_INFO.value, 'test_instance']
+    encoded_request = msgpack.packb(request, use_bin_type=True)
+
+    result = mmp_request_handler.handle_request(encoded_request)
+    decoded_result = msgpack.unpackb(result, raw=False)
+
+    assert decoded_result[0] == ResponseType.SUCCESS.value
+    elapsed_time, checkpoints, resume, snapshot_directory = decoded_result[1:]
+
+    assert elapsed_time > 0.0
+
+    assert isinstance(checkpoints, dict)
+    assert checkpoints.keys() == {'at_end', 'wallclock_time', 'simulation_time'}
+    assert checkpoints['at_end'] is True
+    assert checkpoints['simulation_time'] == []
+    wallclock_time = checkpoints['wallclock_time']
+    assert len(wallclock_time) == 2
+    assert wallclock_time[0] == {'start': None, 'stop': None, 'every': 10}
+    assert wallclock_time[1] == {'at': [1, 2, 3.0]}
+
+    assert resume is not None
+    assert Path(resume) == resume_path
+
+    assert snapshot_directory is None
+
+
+def test_get_checkpoint_info2(registered_mmp_request_handler2, tmp_path):
+    request = [RequestType.GET_CHECKPOINT_INFO.value, 'test_instance']
+    encoded_request = msgpack.packb(request, use_bin_type=True)
+
+    result = registered_mmp_request_handler2.handle_request(encoded_request)
+    decoded_result = msgpack.unpackb(result, raw=False)
+
+    assert decoded_result[0] == ResponseType.SUCCESS.value
+    snapshot_directory = decoded_result[4]
+    assert snapshot_directory == (
+            str(tmp_path) + '/instances/test_instance/snapshots')
+
+
 def test_double_register_instance(mmp_request_handler):
     request = [
             RequestType.REGISTER_INSTANCE.value,
@@ -228,3 +282,24 @@ def test_request_peers_unknown(registered_mmp_request_handler2):
     assert status == ResponseType.ERROR.value
     assert error_msg is not None
     assert 'does_not_exist' in error_msg
+
+
+def test_submit_snapshot(registered_mmp_request_handler):
+    register_snapshot = MagicMock()
+    registered_mmp_request_handler._snapshot_registry.register_snapshot = \
+        register_snapshot
+
+    instance_id = 'micro[1][2]'
+    snapshot = SnapshotMetadata(
+            ['1', '2'], 1.234, 2.345, 3.456,
+            {'in': [1], 'out': [0]}, True, 'fname')
+    snapshot_dict = dataclasses.asdict(snapshot)
+
+    request = [RequestType.SUBMIT_SNAPSHOT.value, instance_id, snapshot_dict]
+    encoded_request = msgpack.packb(request, use_bin_type=True)
+
+    result = registered_mmp_request_handler.handle_request(encoded_request)
+    decoded_result = msgpack.unpackb(result, raw=False)
+
+    assert decoded_result[0] == ResponseType.SUCCESS.value
+    register_snapshot.assert_called_once_with(Reference(instance_id), snapshot)
diff --git a/libmuscle/python/libmuscle/manager/test/test_snapshot_registry.py b/libmuscle/python/libmuscle/manager/test/test_snapshot_registry.py
new file mode 100644
index 00000000..2d71630b
--- /dev/null
+++ b/libmuscle/python/libmuscle/manager/test/test_snapshot_registry.py
@@ -0,0 +1,420 @@
+from datetime import datetime, timedelta
+from pathlib import Path
+from unittest.mock import MagicMock
+
+import pytest
+from ymmsl import (
+        Configuration, Model, Component, Conduit, Implementation,
+        KeepsStateForNextUse, Reference)
+
+from libmuscle.manager.snapshot_registry import (
+    SnapshotNode, SnapshotRegistry, calc_consistency, calc_consistency_list,
+    safe_get, _ConnectionInfo)
+from libmuscle.manager.topology_store import TopologyStore
+from libmuscle.snapshot import SnapshotMetadata
+
+
+def make_snapshot(**msg_counts) -> SnapshotMetadata:
+    return SnapshotMetadata([], 0, 0, 0, {**msg_counts}, False, '')
+
+
+@pytest.fixture(params=[True, False])
+def micro_is_stateless(request: pytest.FixtureRequest) -> bool:
+    return request.param
+
+
+@pytest.fixture
+def macro_micro(micro_is_stateless: bool) -> Configuration:
+    components = [
+            Component('macro', 'macro_impl'),
+            Component('micro', 'micro_impl')]
+    conduits = [
+            Conduit('macro.o_i', 'micro.f_i'),
+            Conduit('micro.o_f', 'macro.s')]
+    model = Model('macro_micro', components, conduits)
+
+    if micro_is_stateless:
+        micro_impl = Implementation(
+                'micro_impl',
+                keeps_state_for_next_use=KeepsStateForNextUse.NO,
+                executable='pass')
+    else:
+        micro_impl = Implementation('micro_impl', executable='pass')
+
+    implementations = [
+            Implementation('macro_impl', executable='pass'),
+            micro_impl]
+
+    return Configuration(model, implementations=implementations)
+
+
+@pytest.fixture
+def uq(macro_micro: Configuration) -> Configuration:
+    for component in macro_micro.model.components:
+        component.multiplicity = [5]
+    macro_micro.model.components.append(Component('qmc', 'qmc_impl'))
+    macro_micro.model.components.append(Component('rr', 'rr_impl'))
+    macro_micro.model.conduits.extend([
+            Conduit('qmc.parameters_out', 'rr.front_in'),
+            Conduit('rr.front_out', 'qmc.states_in'),
+            Conduit('rr.back_out', 'macro.muscle_settings_in'),
+            Conduit('macro.final_state_out', 'rr.back_in')])
+    macro_micro.implementations[Reference('qmc_impl')] = Implementation(
+            'qmc_impl', executable='pass')
+    macro_micro.implementations[Reference('rr_impl')] = Implementation(
+            'rr_impl', executable='pass')
+    return macro_micro
+
+
+def test_safe_get() -> None:
+    assert safe_get([], 0, 1) == 1
+    assert safe_get([3], 0, 1) == 3
+    assert safe_get([3], 1, 5) == 5
+    for i in range(10):
+        expected = -1 if i >= 3 else i + 3
+        assert safe_get([3, 4, 5], i, -1) == expected
+
+
+def test_calc_consistency() -> None:
+    num_sent = 3
+    for num_received in [2, 3, 4, 5]:
+        expect = num_received in [3, 4]
+        assert calc_consistency(num_sent, num_received, True, False) is expect
+        assert calc_consistency(num_received, num_sent, False, False) is expect
+
+    num_received = 10
+    for num_sent in [8, 9, 10, 11]:
+        expect = num_sent in [9, 10]
+        assert calc_consistency(num_sent, num_received, True, False) is expect
+        assert calc_consistency(num_received, num_sent, False, False) is expect
+
+
+def test_calc_consistency_with_restart() -> None:
+    # Check normal rules
+    assert calc_consistency(0, 0, True, True)
+    assert calc_consistency(0, 0, False, True)
+    assert not calc_consistency(1, 0, True, True)
+    assert not calc_consistency(1, 0, True, False)
+    assert calc_consistency(1, 0, False, False)
+    # Different: num2 == 0 comes from the restarted actor, we do not want a
+    # resume file to be created in this instance (because an instance further in
+    # the call chain is ahead of the one that would be restarted):
+    assert not calc_consistency(1, 0, False, True)
+
+
+def test_calc_consistency_list() -> None:
+    num_sent = [3, 3]
+    for num_received in [[2, 3], [3, 2], [3, 5], [], [4, 4, 0, 0, 2]]:
+        assert not calc_consistency_list(num_sent, num_received, True, False)
+        assert not calc_consistency_list(num_received, num_sent, False, False)
+    for num_received in [[3, 3], [3, 4], [4, 3], [4, 4],
+                         [3, 3, 1], [4, 4, 0, 0, 0, 1, 0, 1]]:
+        assert calc_consistency_list(num_sent, num_received, True, False)
+        assert calc_consistency_list(num_received, num_sent, False, False)
+
+
+def test_write_ymmsl(tmp_path: Path):
+    configuration = Configuration(Model('empty', []))
+    snapshot_registry = SnapshotRegistry(
+            configuration, tmp_path, TopologyStore(configuration))
+    snapshot_registry._write_snapshot_ymmsl([])
+
+    paths = list(tmp_path.iterdir())
+    assert len(paths) == 1
+    assert paths[0].suffix == ".ymmsl"
+    paths[0].unlink()
+
+    now = datetime.now()
+    for seconds in range(3):
+        time = (now + timedelta(seconds=seconds)).strftime("%Y%m%d_%H%M%S")
+        (tmp_path / f'snapshot_{time}.ymmsl').touch()
+    snapshot_registry._write_snapshot_ymmsl([])
+    paths = list(tmp_path.iterdir())
+    assert len(paths) == 4
+    paths = list(tmp_path.glob('*_1.ymmsl'))
+    assert len(paths) == 1
+
+
+def test_snapshot_config():
+    configuration = Configuration(Model('empty', []))
+    snapshot_registry = SnapshotRegistry(
+            configuration, None, TopologyStore(configuration))
+    micro_metadata = SnapshotMetadata(
+            ['simulation_time >= 24.0', 'wallclocktime >= 10'],
+            10.123456789, 24.3456789, None, {}, False, 'micro_snapshot')
+    macro_metadata = SnapshotMetadata(
+            ['simulation_time >= 12.0', 'wallclocktime >= 10'],
+            10.123456789, 12.3456789, None, {}, False, 'macro_snapshot')
+    snapshots = [
+            SnapshotNode(1, Reference('micro'), micro_metadata, set()),
+            SnapshotNode(1, Reference('macro'), macro_metadata, set())]
+
+    now = datetime.now()
+    config = snapshot_registry._generate_snapshot_config(snapshots, now)
+    assert len(config.resume) == 2
+    assert config.resume[Reference('macro')] == Path('macro_snapshot')
+    assert config.resume[Reference('micro')] == Path('micro_snapshot')
+    # note: no automatic testing for formatting, should verify by eye if this
+    # looks okay..
+    print(config.description)
+
+    long_metadata = SnapshotMetadata(
+            ['simulation_time >= 24.0'], 1.23456789e-10, 1.23456789e10, None,
+            {}, False, '/this/is/a/long/path/to/the/snapshot/file.pack')
+    snapshots.append(SnapshotNode(
+            1, Reference('this.is.a.long.reference[10]'), long_metadata, set()))
+
+    config = snapshot_registry._generate_snapshot_config(snapshots, now)
+    assert len(config.resume) == 3
+    assert config.resume[Reference('this.is.a.long.reference[10]')] == Path(
+            '/this/is/a/long/path/to/the/snapshot/file.pack')
+    print(config.description)
+
+
+def test_peers(uq: Configuration) -> None:
+    snapshot_registry = SnapshotRegistry(uq, None, TopologyStore(uq))
+    macro = Reference('macro')
+    micro = Reference('micro')
+    qmc = Reference('qmc')
+    rr = Reference('rr')
+
+    all_instances = {qmc, rr} | {macro + i for i in range(5)}
+    all_instances.update(micro + i for i in range(5))
+    assert snapshot_registry._instances == all_instances
+
+    assert snapshot_registry._get_peers(qmc) == {rr}
+    expected_rr_peers = {qmc} | {macro + i for i in range(5)}
+    assert snapshot_registry._get_peers(rr) == expected_rr_peers
+    for i in range(5):
+        assert snapshot_registry._get_peers(macro + i) == {rr, micro + i}
+        assert snapshot_registry._get_peers(micro + i) == {macro + i}
+
+
+def test_connections(uq: Configuration) -> None:
+    snapshot_registry = SnapshotRegistry(uq, None, TopologyStore(uq))
+    macro = Reference('macro')
+    micro = Reference('micro')
+    qmc = Reference('qmc')
+    rr = Reference('rr')
+
+    assert not snapshot_registry._get_connections(qmc, macro + 1)
+    assert not snapshot_registry._get_connections(macro + 3, qmc)
+    assert not snapshot_registry._get_connections(qmc, micro + 0)
+    assert not snapshot_registry._get_connections(micro + 1, qmc)
+    assert not snapshot_registry._get_connections(rr, micro + 4)
+    assert not snapshot_registry._get_connections(micro + 0, rr)
+
+    connections = snapshot_registry._get_connections(rr, qmc)
+    assert len(connections) == 2
+    for rr_port, qmc_port, info in connections:
+        assert rr_port in (Reference('front_out'), Reference('front_in'))
+        assert qmc_port in (Reference('parameters_out'), Reference('states_in'))
+        is_sending = bool(info & _ConnectionInfo.SELF_IS_SENDING)
+        assert is_sending is (rr_port == Reference('front_out'))
+        # Note: actually both are vector ports, but this is undetectable from
+        # the ymmsl configuration. Luckily we treat it the same as scalar-scalar
+        assert not (info & _ConnectionInfo.SELF_IS_VECTOR)
+        assert not (info & _ConnectionInfo.PEER_IS_VECTOR)
+
+    connections = snapshot_registry._get_connections(macro + 0, rr)
+    assert len(connections) == 2
+    for macro_port, rr_port, info in connections:
+        assert macro_port in (
+                Reference('muscle_settings_in'), Reference('final_state_out'))
+        assert rr_port in (Reference('back_out'), Reference('back_in'))
+        is_sending = bool(info & _ConnectionInfo.SELF_IS_SENDING)
+        assert is_sending is (macro_port == Reference('final_state_out'))
+        assert not (info & _ConnectionInfo.SELF_IS_VECTOR)
+        assert (info & _ConnectionInfo.PEER_IS_VECTOR)
+
+    connections = snapshot_registry._get_connections(rr, macro + 1)
+    assert len(connections) == 2
+    for rr_port, macro_port, info in connections:
+        assert macro_port in (
+                Reference('muscle_settings_in'), Reference('final_state_out'))
+        assert rr_port in (Reference('back_out'), Reference('back_in'))
+        is_sending = bool(info & _ConnectionInfo.SELF_IS_SENDING)
+        assert is_sending is (rr_port == Reference('back_out'))
+        assert (info & _ConnectionInfo.SELF_IS_VECTOR)
+        assert not (info & _ConnectionInfo.PEER_IS_VECTOR)
+
+
+def test_implementation(uq: Configuration) -> None:
+    snapshot_registry = SnapshotRegistry(uq, None, TopologyStore(uq))
+
+    qmc_impl = snapshot_registry._implementation(Reference('qmc'))
+    assert qmc_impl.name == 'qmc_impl'
+
+    missing_impl = snapshot_registry._implementation(Reference('missing'))
+    assert missing_impl is None
+
+
+def test_macro_micro_snapshots(macro_micro: Configuration) -> None:
+    snapshot_registry = SnapshotRegistry(
+            macro_micro, None, TopologyStore(macro_micro))
+    # prevent actually writing a ymmsl file, testing that separately
+    snapshot_registry._write_snapshot_ymmsl = MagicMock()
+    macro = Reference('macro')
+    micro = Reference('micro')
+
+    macro_snapshot = make_snapshot(o_i=[3], s=[3])
+    snapshot_registry._add_snapshot(macro, macro_snapshot)
+
+    assert len(snapshot_registry._snapshots[macro]) == 1
+    node = snapshot_registry._snapshots[macro][0]
+    assert node.consistent is False
+    assert node.consistent_peers == {}
+    assert node.instance == macro
+    assert node.num == 1
+    assert node.snapshot is macro_snapshot
+    assert node.peers == {micro}
+    snapshot_registry._write_snapshot_ymmsl.assert_not_called()
+
+    # Note: this snapshot is not realistic, it should have come in before
+    # the macro snapshot above. However, it's still useful for testing the
+    # consistency algorithm
+    micro_snapshot = make_snapshot(f_i=[2], o_f=[1])
+    snapshot_registry._add_snapshot(micro, micro_snapshot)
+
+    assert len(snapshot_registry._snapshots[micro]) == 1
+    assert snapshot_registry._snapshots[micro][0].consistent is False
+    snapshot_registry._write_snapshot_ymmsl.assert_not_called()
+
+    micro_snapshot = make_snapshot(f_i=[3], o_f=[2])
+    snapshot_registry._add_snapshot(micro, micro_snapshot)
+
+    # The first micro snapshots should be cleaned up now
+    assert len(snapshot_registry._snapshots[micro]) == 1
+    micro_node = snapshot_registry._snapshots[micro][0]
+    assert micro_node.consistent
+    snapshot_registry._write_snapshot_ymmsl.assert_called_once_with(
+            [micro_node, node])
+    snapshot_registry._write_snapshot_ymmsl.reset_mock()
+
+    # 3 micro snapshots in the same reuse:
+    for _ in range(3):
+        micro_snapshot = make_snapshot(f_i=[4], o_f=[3])
+        snapshot_registry._add_snapshot(micro, micro_snapshot)
+
+    # Previous micro snapshot should be cleaned up now
+    assert len(snapshot_registry._snapshots[micro]) == 1
+    micro_node = snapshot_registry._snapshots[micro][-1]
+    assert snapshot_registry._write_snapshot_ymmsl.call_count == 3
+    snapshot_registry._write_snapshot_ymmsl.assert_called_with(
+            [micro_node, node])
+    snapshot_registry._write_snapshot_ymmsl.reset_mock()
+
+    macro_snapshot = make_snapshot(o_i=[4], s=[4])
+    snapshot_registry._add_snapshot(macro, macro_snapshot)
+    snapshot_registry._write_snapshot_ymmsl.assert_called_once()
+    snapshot_registry._write_snapshot_ymmsl.reset_mock()
+
+    # 3 micro snapshots in the same reuse, but inconcistent with previous macro
+    for _ in range(3):
+        micro_snapshot = make_snapshot(f_i=[6], o_f=[5])
+        snapshot_registry._add_snapshot(micro, micro_snapshot)
+
+    # All three should be present now in addition to the one last used in
+    # workflow snapshot
+    assert len(snapshot_registry._snapshots[micro]) == 4
+    snapshot_registry._write_snapshot_ymmsl.assert_not_called()
+
+    macro_snapshot = make_snapshot(o_i=[6], s=[6])
+    snapshot_registry._add_snapshot(macro, macro_snapshot)
+    assert snapshot_registry._write_snapshot_ymmsl.call_count == 3
+    assert len(snapshot_registry._snapshots[micro]) == 1
+    assert len(snapshot_registry._snapshots[macro]) == 1
+
+
+def test_uq(uq: Configuration) -> None:
+    snapshot_registry = SnapshotRegistry(uq, None, TopologyStore(uq))
+    # prevent actually writing a ymmsl file, testing that separately
+    snapshot_registry._write_snapshot_ymmsl = MagicMock()
+    macro = Reference('macro')
+    micro = Reference('micro')
+    qmc = Reference('qmc')
+    rr = Reference('rr')
+
+    qmc_snapshot = make_snapshot(parameters_out=[], states_in=[])
+    snapshot_registry._add_snapshot(qmc, qmc_snapshot)
+
+    rr_snapshot = make_snapshot(
+            front_in=[1, 1, 1, 1, 1, 0, 0, 0, 0, 0],
+            front_out=[0] * 10,
+            back_out=[1, 1, 1, 1, 1],
+            back_in=[0] * 5)
+    snapshot_registry._add_snapshot(rr, rr_snapshot)
+    node = snapshot_registry._snapshots[rr][-1]
+    assert qmc in node.consistent_peers
+    snapshot_registry._write_snapshot_ymmsl.assert_not_called()
+
+    macro_snapshot = make_snapshot(
+            muscle_settings_in=[1], final_state_out=[0], o_i=[0], s=[0])
+    for i in range(5):
+        snapshot_registry._add_snapshot(macro + i, macro_snapshot)
+        node = snapshot_registry._snapshots[macro + i][-1]
+        assert node.consistent_peers.keys() == {rr}
+        snapshot_registry._write_snapshot_ymmsl.assert_not_called()
+
+    micro_snapshot = make_snapshot(f_i=[1], o_f=[0])
+    for i in range(5):
+        snapshot_registry._add_snapshot(micro + i, micro_snapshot)
+        node = snapshot_registry._snapshots[micro + i][-1]
+        assert node.consistent_peers.keys() == {macro + i}
+        if i == 4:
+            snapshot_registry._write_snapshot_ymmsl.assert_called_once()
+            snapshot_registry._write_snapshot_ymmsl.reset_mock()
+        else:
+            snapshot_registry._write_snapshot_ymmsl.assert_not_called()
+
+    qmc_snapshot = make_snapshot(parameters_out=[1, 1, 1, 1, 1], states_in=[])
+    snapshot_registry._add_snapshot(qmc, qmc_snapshot)
+    node = snapshot_registry._snapshots[qmc][-1]
+    assert node.consistent_peers.keys() == {rr}
+    snapshot_registry._write_snapshot_ymmsl.assert_called_once()
+    snapshot_registry._write_snapshot_ymmsl.reset_mock()
+    assert len(snapshot_registry._snapshots[qmc]) == 1  # previous is cleaned up
+
+
+def test_heuristic_rollbacks() -> None:
+    components = [Component(f'comp{i}', f'impl{i}') for i in range(4)]
+    conduits = [Conduit(f'comp{i}.o_f', f'comp{i+1}.f_i') for i in range(3)]
+    model = Model('linear', components, conduits)
+    implementations = [
+            Implementation(f'impl{i}', script='xyz')
+            for i in range(4)]
+    config = Configuration(model, implementations=implementations)
+
+    comp1, comp2, comp3, comp4 = (Reference(f'comp{i}') for i in range(4))
+
+    snapshot_registry = SnapshotRegistry(config, None, TopologyStore(config))
+    # prevent actually writing a ymmsl file, testing that separately
+    snapshot_registry._write_snapshot_ymmsl = MagicMock()
+
+    for i in range(4):
+        snapshot_registry._add_snapshot(comp1, make_snapshot(o_f=[i]))
+    assert len(snapshot_registry._snapshots[comp1]) == 4
+
+    for i in range(10):
+        snapshot_registry._add_snapshot(
+                comp2, make_snapshot(f_i=[1], o_f=[0]))
+        snapshot_registry._add_snapshot(
+                comp3, make_snapshot(f_i=[1], o_f=[0]))
+    assert len(snapshot_registry._snapshots[comp2]) == 10
+    assert len(snapshot_registry._snapshots[comp3]) == 10
+
+    snapshot_registry._add_snapshot(comp2, make_snapshot(f_i=[2], o_f=[1]))
+    assert len(snapshot_registry._snapshots[comp2]) == 11
+    snapshot_registry._add_snapshot(comp2, make_snapshot(f_i=[3], o_f=[2]))
+    assert len(snapshot_registry._snapshots[comp2]) == 12
+
+    snapshot_registry._write_snapshot_ymmsl.assert_not_called()
+
+    snapshot_registry._add_snapshot(
+            comp4, make_snapshot(f_i=[1]))
+    snapshot_registry._write_snapshot_ymmsl.assert_called()
+
+    assert len(snapshot_registry._snapshots[comp1]) == 2
+    assert len(snapshot_registry._snapshots[comp2]) == 2
+    assert len(snapshot_registry._snapshots[comp3]) == 1
+    assert len(snapshot_registry._snapshots[comp4]) == 1
diff --git a/libmuscle/python/libmuscle/test/conftest.py b/libmuscle/python/libmuscle/test/conftest.py
index a8d0ad72..77422ee9 100644
--- a/libmuscle/python/libmuscle/test/conftest.py
+++ b/libmuscle/python/libmuscle/test/conftest.py
@@ -3,6 +3,7 @@
 
 from ymmsl import Settings
 
+from libmuscle.api_guard import APIGuard
 from libmuscle.communicator import Message
 from libmuscle.mmp_client import MMPClient
 
@@ -21,3 +22,8 @@ def message() -> Message:
 @pytest.fixture
 def message2() -> Message:
     return Message(0.0, None, {'test': 17}, Settings())
+
+
+@pytest.fixture
+def guard() -> APIGuard:
+    return APIGuard()
diff --git a/libmuscle/python/libmuscle/test/test_api_guard.py b/libmuscle/python/libmuscle/test/test_api_guard.py
new file mode 100644
index 00000000..f67bde93
--- /dev/null
+++ b/libmuscle/python/libmuscle/test/test_api_guard.py
@@ -0,0 +1,154 @@
+from typing import Callable, Set
+
+import pytest
+
+from libmuscle.api_guard import APIGuard
+
+
+def test_no_checkpointing_support(guard):
+    for _ in range(3):
+        guard.verify_reuse_instance()
+        guard.reuse_instance_done(True)
+
+    assert not guard.uses_checkpointing()
+
+    guard.verify_reuse_instance()
+    guard.reuse_instance_done(False)
+
+    assert not guard.uses_checkpointing()
+
+
+def test_final_snapshot_only(guard):
+    for i in range(4):
+        guard.verify_reuse_instance()
+        guard.reuse_instance_done(True)
+
+        guard.verify_resuming()
+        if i == 0:
+            guard.resuming_done(True)
+
+            guard.verify_load_snapshot()
+            guard.load_snapshot_done()
+        else:
+            guard.resuming_done(False)
+
+        guard.verify_should_init()
+        guard.should_init_done()
+
+        guard.verify_should_save_final_snapshot()
+        if i == 2:
+            guard.should_save_final_snapshot_done(True)
+
+            guard.verify_save_final_snapshot()
+            guard.save_final_snapshot_done()
+        else:
+            guard.should_save_final_snapshot_done(False)
+
+    guard.verify_reuse_instance()
+    guard.reuse_instance_done(False)
+
+
+def test_full_checkpointing(guard):
+    for i in range(4):
+        guard.verify_reuse_instance()
+        guard.reuse_instance_done(True)
+
+        guard.verify_resuming()
+        if i == 0:
+            guard.resuming_done(True)
+
+            guard.verify_load_snapshot()
+            guard.load_snapshot_done()
+        else:
+            guard.resuming_done(False)
+
+        guard.verify_should_init()
+        guard.should_init_done()
+
+        for j in range(3):
+            guard.verify_should_save_snapshot()
+            if j != 2:
+                guard.should_save_snapshot_done(True)
+
+                guard.verify_save_snapshot()
+                guard.save_snapshot_done()
+            else:
+                guard.should_save_snapshot_done(False)
+
+        guard.verify_should_save_final_snapshot()
+        if i == 2:
+            guard.should_save_final_snapshot_done(True)
+
+            guard.verify_save_final_snapshot()
+            guard.save_final_snapshot_done()
+        else:
+            guard.should_save_final_snapshot_done(False)
+
+    guard.verify_reuse_instance()
+    guard.reuse_instance_done(False)
+
+
+_api_guard_funs = (
+    (APIGuard.verify_reuse_instance, ()),
+    (APIGuard.reuse_instance_done, (True,)),
+    (APIGuard.verify_resuming, ()),
+    (APIGuard.resuming_done, (True,)),
+    (APIGuard.verify_load_snapshot, ()),
+    (APIGuard.load_snapshot_done, ()),
+    (APIGuard.verify_should_init, ()),
+    (APIGuard.should_init_done, ()),
+    (APIGuard.verify_should_save_snapshot, ()),
+    (APIGuard.should_save_snapshot_done, (True,)),
+    (APIGuard.verify_save_snapshot, ()),
+    (APIGuard.save_snapshot_done, ()),
+    (APIGuard.verify_should_save_final_snapshot, ()),
+    (APIGuard.should_save_final_snapshot_done, (True,)),
+    (APIGuard.verify_save_final_snapshot, ())
+)
+
+
+def run_until_before(guard: APIGuard, excluded: Callable) -> None:
+    for fun, args in _api_guard_funs:
+        if fun is excluded:
+            break
+        fun(guard, *args)
+
+
+def check_all_raise_except(guard: APIGuard, excluded: Set[Callable]) -> None:
+    for fun, args in _api_guard_funs:
+        if fun.__name__.startswith('verify_'):
+            if fun not in excluded:
+                with pytest.raises(RuntimeError):
+                    fun(guard, *args)
+            else:
+                fun(guard, *args)
+
+
+@pytest.mark.parametrize('fun', [
+        APIGuard.verify_load_snapshot,
+        APIGuard.verify_should_init, APIGuard.verify_save_snapshot,
+        APIGuard.verify_save_final_snapshot])
+def test_missing_step(guard, fun):
+    run_until_before(guard, fun)
+    check_all_raise_except(guard, {fun})
+
+
+def test_missing_resuming(guard):
+    run_until_before(guard, APIGuard.verify_resuming)
+    check_all_raise_except(guard, {
+        APIGuard.verify_resuming, APIGuard.verify_reuse_instance})
+
+
+def test_missing_should_save_final(guard):
+    run_until_before(guard, APIGuard.verify_should_save_final_snapshot)
+    check_all_raise_except(guard, {
+        APIGuard.verify_should_save_snapshot,
+        APIGuard.verify_should_save_final_snapshot})
+
+
+def test_double_should_save(guard):
+    run_until_before(guard, APIGuard.verify_should_save_snapshot)
+    guard.verify_should_save_snapshot()
+    guard.should_save_snapshot_done(True)
+    with pytest.raises(RuntimeError):
+        guard.verify_should_save_snapshot()
diff --git a/libmuscle/python/libmuscle/test/test_checkpoint_triggers.py b/libmuscle/python/libmuscle/test/test_checkpoint_triggers.py
new file mode 100644
index 00000000..e111a758
--- /dev/null
+++ b/libmuscle/python/libmuscle/test/test_checkpoint_triggers.py
@@ -0,0 +1,186 @@
+import time
+import pytest
+from ymmsl import CheckpointRangeRule, CheckpointAtRule, Checkpoints
+
+from libmuscle.checkpoint_triggers import (
+    CombinedCheckpointTriggers, AtCheckpointTrigger, RangeCheckpointTrigger,
+    TriggerManager)
+
+
+def test_at_checkpoint_trigger():
+    trigger = AtCheckpointTrigger([CheckpointAtRule([1, 3, 4, 4.5, 9])])
+
+    assert trigger.next_checkpoint(0) == 1
+    assert trigger.previous_checkpoint(0) is None
+
+    assert trigger.next_checkpoint(1) == 3
+    assert trigger.previous_checkpoint(1) == 1
+
+    eps = 1e-16
+    assert trigger.next_checkpoint(1 - eps) == 1
+    assert trigger.previous_checkpoint(1 - eps) is None
+
+    assert trigger.next_checkpoint(3.9) == 4
+    assert trigger.previous_checkpoint(3.9) == 3
+
+    assert trigger.next_checkpoint(4.1) == 4.5
+    assert trigger.previous_checkpoint(4.1) == 4
+
+    assert trigger.next_checkpoint(5) == 9
+    assert trigger.previous_checkpoint(5) == 4.5
+
+    assert trigger.next_checkpoint(9) is None
+    assert trigger.previous_checkpoint(9) == 9
+
+    assert trigger.next_checkpoint(11) is None
+    assert trigger.previous_checkpoint(11) == 9
+
+
+def test_range_checkpoint_trigger():
+    range = CheckpointRangeRule(start=0, stop=20, every=1.2)
+    trigger = RangeCheckpointTrigger(range)
+
+    assert trigger.next_checkpoint(-1) == 0
+    assert trigger.previous_checkpoint(-1) is None
+
+    assert trigger.next_checkpoint(0) == pytest.approx(1.2)
+    assert trigger.previous_checkpoint(0) == 0
+
+    assert trigger.next_checkpoint(8) == pytest.approx(8.4)
+    assert trigger.previous_checkpoint(8) == pytest.approx(7.2)
+
+    assert trigger.next_checkpoint(18.2) == pytest.approx(19.2)
+    assert trigger.previous_checkpoint(18.2) == pytest.approx(18)
+
+    assert trigger.next_checkpoint(20) is None
+    assert trigger.previous_checkpoint(20) == pytest.approx(19.2)
+
+
+def test_range_checkpoint_trigger_default_stop():
+    range = CheckpointRangeRule(start=1, every=1.2)
+    trigger = RangeCheckpointTrigger(range)
+
+    assert trigger.next_checkpoint(-1.) == 1
+    assert trigger.previous_checkpoint(-1.) is None
+
+    assert trigger.next_checkpoint(148148.) == pytest.approx(148148.2)
+    assert trigger.previous_checkpoint(148148.) == pytest.approx(148147)
+
+    assert trigger.next_checkpoint(148148148.) == pytest.approx(148148149)
+    assert trigger.previous_checkpoint(148148148.) == pytest.approx(148148147.8)
+
+
+def test_range_checkpoint_trigger_default_start():
+    range = CheckpointRangeRule(every=1.2, stop=10)
+    trigger = RangeCheckpointTrigger(range)
+
+    assert trigger.next_checkpoint(10) is None
+    assert trigger.previous_checkpoint(10) == pytest.approx(9.6)
+
+    assert trigger.next_checkpoint(0.0) == pytest.approx(1.2)
+    assert trigger.previous_checkpoint(0.0) == pytest.approx(0.0)
+
+    assert trigger.next_checkpoint(-148148.) == pytest.approx(-148147.2)
+    assert trigger.previous_checkpoint(-148148.) == pytest.approx(-148148.4)
+
+
+def test_combined_checkpoint_trigger_every_at():
+    rules = [CheckpointRangeRule(every=10), CheckpointAtRule([3, 7, 13, 17])]
+    trigger = CombinedCheckpointTriggers(rules)
+
+    assert trigger.next_checkpoint(-11.) == pytest.approx(-10)
+    assert trigger.previous_checkpoint(-11) == pytest.approx(-20)
+
+    assert trigger.next_checkpoint(0.) == pytest.approx(3)
+    assert trigger.previous_checkpoint(0.) == pytest.approx(0)
+
+    assert trigger.next_checkpoint(8.3) == pytest.approx(10)
+    assert trigger.previous_checkpoint(8.3) == pytest.approx(7)
+
+    assert trigger.next_checkpoint(14.2) == pytest.approx(17)
+    assert trigger.previous_checkpoint(14.2) == pytest.approx(13)
+
+    assert trigger.next_checkpoint(25.2) == pytest.approx(30)
+    assert trigger.previous_checkpoint(25.2) == pytest.approx(20)
+
+
+def test_combined_checkpoint_trigger_at_ranges():
+    rules = [CheckpointAtRule([3, 7, 13, 17]),
+             CheckpointRangeRule(start=0, every=5, stop=20),
+             CheckpointRangeRule(start=20, every=20, stop=100)]
+    trigger = CombinedCheckpointTriggers(rules)
+
+    assert trigger.next_checkpoint(-11.) == pytest.approx(0)
+    assert trigger.previous_checkpoint(-11) is None
+
+    assert trigger.next_checkpoint(0.) == pytest.approx(3)
+    assert trigger.previous_checkpoint(0.) == pytest.approx(0)
+
+    assert trigger.next_checkpoint(8.3) == pytest.approx(10)
+    assert trigger.previous_checkpoint(8.3) == pytest.approx(7)
+
+    assert trigger.next_checkpoint(14.2) == pytest.approx(15)
+    assert trigger.previous_checkpoint(14.2) == pytest.approx(13)
+
+    assert trigger.next_checkpoint(19.3) == pytest.approx(20)
+    assert trigger.previous_checkpoint(19.3) == pytest.approx(17)
+
+    assert trigger.next_checkpoint(25.2) == pytest.approx(40)
+    assert trigger.previous_checkpoint(25.2) == pytest.approx(20)
+
+    assert trigger.next_checkpoint(95.2) == pytest.approx(100)
+    assert trigger.previous_checkpoint(95.2) == pytest.approx(80)
+
+    assert trigger.next_checkpoint(125.2) is None
+    assert trigger.previous_checkpoint(125.2) == pytest.approx(100)
+
+
+def test_trigger_manager_reference_time():
+    monotonic_start = time.monotonic()
+    ref_elapsed = 15.0
+    trigger_manager = TriggerManager()
+    trigger_manager.set_checkpoint_info(ref_elapsed, Checkpoints(at_end=True))
+    elapsed_walltime = trigger_manager.elapsed_walltime()
+    duration = time.monotonic() - monotonic_start
+    assert ref_elapsed < elapsed_walltime <= (ref_elapsed + duration)
+
+
+def test_trigger_manager():
+    ref_elapsed = 0.0
+    trigger_manager = TriggerManager()
+    trigger_manager.set_checkpoint_info(ref_elapsed, Checkpoints(
+            at_end=True,
+            wallclock_time=[CheckpointAtRule([1e-12])],
+            simulation_time=[CheckpointAtRule([1, 3, 5])]))
+
+    assert trigger_manager.should_save_snapshot(0.1)
+    triggers = trigger_manager.get_triggers()
+    assert len(triggers) == 1
+    assert "wallclock_time" in triggers[0]
+    trigger_manager.update_checkpoints(0.1)
+
+    assert not trigger_manager.should_save_snapshot(0.99)
+
+    assert trigger_manager.should_save_snapshot(3.2)
+    triggers = trigger_manager.get_triggers()
+    assert len(triggers) == 1
+    assert "simulation_time" in triggers[0]
+    trigger_manager.update_checkpoints(3.2)
+
+    assert trigger_manager.should_save_final_snapshot(True, 7.0)
+    assert len(trigger_manager.get_triggers()) > 0
+    trigger_manager.update_checkpoints(7.0)
+
+    assert not trigger_manager.should_save_snapshot(7.1)
+
+    assert trigger_manager.should_save_final_snapshot(False, None)
+    trigger_manager.update_checkpoints(7.1)
+
+
+def test_no_checkpointing() -> None:
+    trigger_manager = TriggerManager()
+    trigger_manager.set_checkpoint_info(0.0, Checkpoints())
+
+    assert not trigger_manager.should_save_snapshot(1)
+    assert not trigger_manager.should_save_snapshot(5000)
+    assert not trigger_manager.should_save_final_snapshot(False, None)
diff --git a/libmuscle/python/libmuscle/test/test_communicator.py b/libmuscle/python/libmuscle/test/test_communicator.py
index e2414856..140e9399 100644
--- a/libmuscle/python/libmuscle/test/test_communicator.py
+++ b/libmuscle/python/libmuscle/test/test_communicator.py
@@ -1,3 +1,4 @@
+import logging
 from typing import List
 from libmuscle.communicator import Communicator, Endpoint, Message
 from libmuscle.mpp_message import ClosePort, MPPMessage
@@ -68,6 +69,8 @@ def gpe(p, s) -> List[Reference]:
     communicator._ports = {
             'out': Port('out', Operator.O_I, False, True, 1, []),
             'in': Port('in', Operator.S, False, True, 1, [])}
+    communicator._muscle_settings_in = \
+        communicator._Communicator__settings_in_port([])
     yield communicator
     communicator.shutdown()
 
@@ -104,6 +107,8 @@ def gpe(p, s) -> List[Reference]:
     communicator._ports = {
             'out': Port('out', Operator.O_I, True, True, 0, [20]),
             'in': Port('in', Operator.S, True, True, 0, [20])}
+    communicator._muscle_settings_in = \
+        communicator._Communicator__settings_in_port([])
     yield communicator
     communicator.shutdown()
 
@@ -282,6 +287,7 @@ def test_send_message(communicator, message) -> None:
     assert msg.timestamp == 0.0
     assert msg.next_timestamp is None
     assert msg.settings_overlay == Settings()
+    assert msg.message_number == 0
     assert msg.data == b'test'
 
 
@@ -305,6 +311,7 @@ def test_send_msgpack(communicator, message2) -> None:
     assert msg.sender == 'kernel[13].out'
     assert msg.receiver == 'other.in[13]'
     assert msg.settings_overlay == Settings()
+    assert msg.message_number == 0
     assert msg.data == {'test': 17}
 
 
@@ -319,6 +326,7 @@ def test_send_message_with_slot(communicator2, message) -> None:
     assert msg.sender == 'other.out[13]'
     assert msg.receiver == 'kernel[13].in'
     assert msg.settings_overlay == Settings()
+    assert msg.message_number == 0
     assert msg.data == b'test'
 
 
@@ -349,6 +357,7 @@ def test_send_message_with_settings(communicator, message) -> None:
     assert msg.sender == 'kernel[13].out'
     assert msg.receiver == 'other.in[13]'
     assert msg.settings_overlay.as_ordered_dict() == {'test2': 'testing'}
+    assert msg.message_number == 0
     assert msg.data == b'test'
 
 
@@ -364,6 +373,7 @@ def test_send_settings(communicator, message) -> None:
     assert msg.sender == 'kernel[13].out'
     assert msg.receiver == 'other.in[13]'
     assert msg.settings_overlay == Settings()
+    assert msg.message_number == 0
     assert msg.data == Settings({'test1': 'testing'})
 
 
@@ -379,6 +389,7 @@ def test_close_port(communicator) -> None:
     assert msg.timestamp == float('inf')
     assert msg.next_timestamp is None
     assert msg.settings_overlay == Settings()
+    assert msg.message_number == 0
     assert isinstance(msg.data, ClosePort)
 
 
@@ -386,28 +397,31 @@ def test_receive_message(communicator) -> None:
     client_mock = MagicMock()
     client_mock.receive.return_value = MPPMessage(
             Reference('other.out[13]'), Reference('kernel[13].in'),
-            None, 0.0, None, Settings({'test1': 12}),
+            None, 0.0, None, Settings({'test1': 12}), 0, 2.0,
             b'test').encoded()
     get_client_mock = MagicMock(return_value=client_mock)
     communicator._Communicator__get_client = get_client_mock
     communicator._profiler = MagicMock()
 
-    msg = communicator.receive_message('in')
+    msg, last_saved = communicator.receive_message('in')
 
     get_client_mock.assert_called_with(Reference('other'))
     client_mock.receive.assert_called_with(Reference('kernel[13].in'))
     assert msg.data == b'test'
     assert msg.settings['test1'] == 12
+    assert last_saved == 2.0
 
 
 def test_receive_message_default(communicator) -> None:
     communicator._peer_manager.is_connected.return_value = False
     default_msg = Message(3.0, 4.0, 'test', Settings())
-    msg = communicator.receive_message('not_connected', default=default_msg)
+    msg, last_saved = communicator.receive_message(
+            'not_connected', default=default_msg)
     assert msg.timestamp == 3.0
     assert msg.next_timestamp == 4.0
     assert msg.data == 'test'
     assert len(msg.settings) == 0
+    assert last_saved == float('-inf')
 
 
 def test_receive_message_no_default(communicator) -> None:
@@ -425,71 +439,75 @@ def test_receive_msgpack(communicator) -> None:
     client_mock = MagicMock()
     client_mock.receive.return_value = MPPMessage(
             Reference('other.out[13]'), Reference('kernel[13].in'),
-            None, 0.0, None, Settings({'test1': 12}),
+            None, 0.0, None, Settings({'test1': 12}), 0, 1.0,
             {'test': 13}).encoded()
     get_client_mock = MagicMock(return_value=client_mock)
     communicator._Communicator__get_client = get_client_mock
     communicator._profiler = MagicMock()
 
-    msg = communicator.receive_message('in')
+    msg, last_saved = communicator.receive_message('in')
 
     get_client_mock.assert_called_with(Reference('other'))
     client_mock.receive.assert_called_with(Reference('kernel[13].in'))
     assert msg.data == {'test': 13}
+    assert last_saved == 1.0
 
 
 def test_receive_with_slot(communicator2) -> None:
     client_mock = MagicMock()
     client_mock.receive.return_value = MPPMessage(
             Reference('kernel[13].out'), Reference('other.in[13]'),
-            None, 0.0, None, Settings({'test': 'testing'}),
+            None, 0.0, None, Settings({'test': 'testing'}), 0, 3.0,
             b'test').encoded()
     get_client_mock = MagicMock(return_value=client_mock)
     communicator2._Communicator__get_client = get_client_mock
     communicator2._profiler = MagicMock()
 
-    msg = communicator2.receive_message('in', 13)
+    msg, last_saved = communicator2.receive_message('in', 13)
 
     get_client_mock.assert_called_with(Reference('kernel[13]'))
     client_mock.receive.assert_called_with(Reference('other.in[13]'))
     assert msg.data == b'test'
     assert msg.settings['test'] == 'testing'
+    assert last_saved == 3.0
 
 
 def test_receive_message_resizable(communicator3) -> None:
     client_mock = MagicMock()
     client_mock.receive.return_value = MPPMessage(
             Reference('other.out[13]'), Reference('kernel.in[13]'),
-            20, 0.0, None, Settings({'test': 'testing'}),
+            20, 0.0, None, Settings({'test': 'testing'}), 0, 12.3,
             b'test').encoded()
     get_client_mock = MagicMock(return_value=client_mock)
     communicator3._Communicator__get_client = get_client_mock
     communicator3._profiler = MagicMock()
 
-    msg = communicator3.receive_message('in', 13)
+    msg, last_saved = communicator3.receive_message('in', 13)
 
     get_client_mock.assert_called_with(Reference('other'))
     client_mock.receive.assert_called_with(Reference('kernel.in[13]'))
     assert msg.data == b'test'
     assert communicator3.get_port('in').get_length() == 20
+    assert last_saved == 12.3
 
 
 def test_receive_with_settings(communicator) -> None:
     client_mock = MagicMock()
     client_mock.receive.return_value = MPPMessage(
             Reference('other.out[13]'), Reference('kernel[13].in'),
-            None, 0.0, None, Settings({'test2': 3.1}),
+            None, 0.0, None, Settings({'test2': 3.1}), 0, 0.1,
             b'test').encoded()
     get_client_mock = MagicMock(return_value=client_mock)
     communicator._Communicator__get_client = get_client_mock
     communicator._profiler = MagicMock()
 
-    msg = communicator.receive_message('in')
+    msg, last_saved = communicator.receive_message('in')
 
     get_client_mock.assert_called_with(Reference('other'))
     client_mock.receive.assert_called_with(Reference('kernel[13].in'))
     assert msg.data == b'test'
     assert msg.settings['test2'] == 3.1
+    assert last_saved == 0.1
 
 
 def test_receive_msgpack_with_slot_and_settings(communicator2) -> None:
@@ -497,55 +515,179 @@ def test_receive_msgpack_with_slot_and_settings(communicator2) -> None:
     client_mock.receive.return_value = MPPMessage(
             Reference('kernel[13].out'), Reference('other.in[13]'),
             None, 0.0, 1.0,
-            Settings({'test': 'testing'}), 'test').encoded()
+            Settings({'test': 'testing'}), 0, 1.0, 'test').encoded()
     get_client_mock = MagicMock(return_value=client_mock)
     communicator2._Communicator__get_client = get_client_mock
     communicator2._profiler = MagicMock()
 
-    msg = communicator2.receive_message('in', 13)
+    msg, last_saved = communicator2.receive_message('in', 13)
 
     get_client_mock.assert_called_with(Reference('kernel[13]'))
     client_mock.receive.assert_called_with(Reference('other.in[13]'))
     assert msg.data == 'test'
     assert msg.settings['test'] == 'testing'
+    assert last_saved == 1.0
 
 
 def test_receive_settings(communicator) -> None:
     client_mock = MagicMock()
     client_mock.receive.return_value = MPPMessage(
             Reference('other.out[13]'), Reference('kernel[13].in'),
-            None, 0.0, None, Settings({'test1': 12}),
+            None, 0.0, None, Settings({'test1': 12}), 0, 1.0,
             Settings({'test': 13})).encoded()
     get_client_mock = MagicMock(return_value=client_mock)
     communicator._Communicator__get_client = get_client_mock
     communicator._profiler = MagicMock()
 
-    msg = communicator.receive_message('in')
+    msg, last_saved = communicator.receive_message('in')
 
     get_client_mock.assert_called_with(Reference('other'))
     client_mock.receive.assert_called_with(Reference('kernel[13].in'))
     assert isinstance(msg.data, Settings)
     assert msg.data['test'] == 13
+    assert last_saved == 1.0
 
 
 def test_receive_close_port(communicator) -> None:
     client_mock = MagicMock()
     client_mock.receive.return_value = MPPMessage(
             Reference('other.out[13]'), Reference('kernel[13].in'),
-            None, 0.0, None, Settings(), ClosePort()).encoded()
+            None, 0.0, None, Settings(), 0, 0.1, ClosePort()).encoded()
     get_client_mock = MagicMock(return_value=client_mock)
     communicator._Communicator__get_client = get_client_mock
     communicator._profiler = MagicMock()
 
-    msg = communicator.receive_message('in')
+    msg, _ = communicator.receive_message('in')
 
     assert isinstance(msg.data, ClosePort)
 
 
 def test_get_message(communicator, message) -> None:
-    communicator.send_message('out', message)
+    communicator.send_message('out', message, None, 2.0)
     ref_message = MPPMessage(
             Reference('kernel[13].out'), Reference('other.in[13]'),
-            None, 0.0, None, Settings(), b'test').encoded()
+            None, 0.0, None, Settings(), 0, 2.0, b'test').encoded()
     assert communicator._post_office.get_message(
             'other.in[13]') == ref_message
+
+
+def test_port_message_counts(communicator, message) -> None:
+    communicator.send_message('out', message)
+    msg_counts = communicator.get_message_counts()
+    assert msg_counts == {'out': [1],
+                          'in': [0],
+                          'muscle_settings_in': [0]}
+
+    communicator.restore_message_counts({'out': [3],
+                                         'in': [2],
+                                         'muscle_settings_in': [4]})
+    communicator.send_message('out', message)
+    msg_counts = communicator.get_message_counts()
+    assert msg_counts == {'out': [4],
+                          'in': [2],
+                          'muscle_settings_in': [4]}
+
+    # empty post office
+    communicator._post_office.get_message('other.in[13]')
+    communicator._post_office.get_message('other.in[13]')
+
+    with pytest.raises(RuntimeError):
+        communicator.restore_message_counts({"x?invalid_port": 3})
+
+
+def test_vector_port_message_counts(communicator2, message) -> None:
+    msg_counts = communicator2.get_message_counts()
+    assert msg_counts == {'out': [0] * 20,
+                          'in': [0] * 20,
+                          'muscle_settings_in': [0]}
+
+    communicator2.send_message('out', message, 13)
+    msg_counts = communicator2.get_message_counts()
+    assert msg_counts == {'out': [0] * 13 + [1] + [0] * 6,
+                          'in': [0] * 20,
+                          'muscle_settings_in': [0]}
+
+    communicator2.restore_message_counts({'out': list(range(20)),
+                                          'in': list(range(20)),
+                                          'muscle_settings_in': [4]})
+    communicator2.send_message('out', message, 13)
+    msg_counts = communicator2.get_message_counts()
+    assert msg_counts == {'out': list(range(13)) + [14] + list(range(14, 20)),
+                          'in': list(range(20)),
+                          'muscle_settings_in': [4]}
+
+    # empty post office
+    communicator2._post_office.get_message('kernel[13].in')
+    communicator2._post_office.get_message('kernel[13].in')
+
+
+def test_port_count_validation(communicator):
+    client_mock = MagicMock()
+    client_mock.receive.return_value = MPPMessage(
+            Reference('other.out[13]'), Reference('kernel[13].in'),
+            None, 0.0, None, Settings({'test1': 12}), 0, 7.6,
+            b'test').encoded()
+    get_client_mock = MagicMock(return_value=client_mock)
+    communicator._Communicator__get_client = get_client_mock
+    communicator._profiler = MagicMock()
+
+    communicator.receive_message('in')
+    assert communicator.get_message_counts()['in'] == [1]
+
+    with pytest.raises(RuntimeError):
+        # the message received has message_number = 0 again
+        communicator.receive_message('in')
+
+
+def test_port_discard_error_on_resume(caplog, communicator):
+    client_mock = MagicMock()
+    client_mock.receive.return_value = MPPMessage(
+            Reference('other.out[13]'), Reference('kernel[13].in'),
+            None, 0.0, None, Settings({'test1': 12}), 1, 2.3,
+            b'test').encoded()
+    get_client_mock = MagicMock(return_value=client_mock)
+    communicator._Communicator__get_client = get_client_mock
+    communicator._profiler = MagicMock()
+
+    communicator.restore_message_counts({'out': [0],
+                                         'in': [2],
+                                         'muscle_settings_in': [0]})
+    for port in communicator._ports.values():
+        assert port._is_resuming == [True]
+        assert port.is_resuming(None)
+
+    # In the next block, the first message with message_number=1 is discarded.
+    # The RuntimeError is raised when 'receiving' the second message with
+    # message_number=1
+    with caplog.at_level(logging.DEBUG, 'libmuscle.communicator'):
+        with pytest.raises(RuntimeError):
+            communicator.receive_message('in')
+        # records 0, 2 and 3 are debug logs for starting/receiving on port
+        assert 'Discarding received message' in caplog.records[1].message
+
+
+def test_port_discard_success_on_resume(caplog, communicator):
+    client_mock = MagicMock()
+    client_mock.receive.side_effect = [MPPMessage(
+            Reference('other.out[13]'), Reference('kernel[13].in'),
+            None, 0.0, None, Settings({'test1': 12}), message_number, 1.0,
+            {'this is message': message_number}).encoded()
+            for message_number in [1, 2]]
+    get_client_mock = MagicMock(return_value=client_mock)
+    communicator._Communicator__get_client = get_client_mock
+    communicator._profiler = MagicMock()
+
+    communicator.restore_message_counts({'out': [0],
+                                         'in': [2],
+                                         'muscle_settings_in': [0]})
+    for port in communicator._ports.values():
+        assert port._is_resuming == [True]
+        assert port.is_resuming(None)
+
+    with caplog.at_level(logging.DEBUG, 'libmuscle.communicator'):
+        msg, _ = communicator.receive_message('in')
+        # records 0, 2 and 3 are debug logs for starting/receiving on port
+        assert 'Discarding received message' in caplog.records[1].message
+    # message_number=1 should be discarded:
+    assert msg.data == {'this is message': 2}
+    assert communicator.get_message_counts()['in'] == [3]
diff --git a/libmuscle/python/libmuscle/test/test_instance.py b/libmuscle/python/libmuscle/test/test_instance.py
index 99d8d37b..7c8b1be4 100644
--- a/libmuscle/python/libmuscle/test/test_instance.py
+++ b/libmuscle/python/libmuscle/test/test_instance.py
@@ -3,7 +3,7 @@
 from unittest.mock import MagicMock, patch
 
 import pytest
-from ymmsl import Operator, Reference, Settings
+from ymmsl import Operator, Reference, Settings, Checkpoints
 
 from libmuscle.communicator import Message
 from libmuscle.instance import Instance
@@ -36,18 +36,20 @@ def sys_argv_instance() -> Generator[None, None, None]:
 
 
 @pytest.fixture
-def instance(sys_argv_instance):
+def instance(sys_argv_instance, tmp_path):
     with patch('libmuscle.instance.MMPClient') as mmp_client, \
          patch('libmuscle.instance.Communicator') as comm_type:
         communicator = MagicMock()
         settings = Settings()
         settings['test1'] = 12
         msg = Message(0.0, 1.0, 'message', settings)
-        communicator.receive_message.return_value = msg
+        communicator.receive_message.return_value = msg, 10.0
         comm_type.return_value = communicator
 
         mmp_client_object = MagicMock()
         mmp_client_object.request_peers.return_value = (None, None, None)
+        checkpoint_info = (0.0, Checkpoints(), None, tmp_path)
+        mmp_client_object.get_checkpoint_info.return_value = checkpoint_info
         mmp_client.return_value = mmp_client_object
 
         instance = Instance({
@@ -59,11 +61,13 @@ def instance(sys_argv_instance):
 
 
 @pytest.fixture
-def instance2(sys_argv_instance):
+def instance2(sys_argv_instance, tmp_path):
     with patch('libmuscle.instance.MMPClient') as mmp_client, \
          patch('libmuscle.instance.Communicator'):
         mmp_client_object = MagicMock()
         mmp_client_object.request_peers.return_value = (None, None, None)
+        checkpoint_info = (0.0, Checkpoints(), None, tmp_path)
+        mmp_client_object.get_checkpoint_info.return_value = checkpoint_info
         mmp_client.return_value = mmp_client_object
         instance = Instance({
             Operator.F_INIT: ['in[]'],
@@ -72,11 +76,13 @@ def instance2(sys_argv_instance):
 
 
 def test_create_instance(
-        sys_argv_instance, log_file_in_tmpdir, sys_argv_manager):
+        sys_argv_instance, log_file_in_tmpdir, sys_argv_manager, tmp_path):
     with patch('libmuscle.instance.MMPClient') as mmp_client, \
          patch('libmuscle.instance.Communicator') as comm_type:
         mmp_client_object = MagicMock()
         mmp_client_object.request_peers.return_value = (None, None, None)
+        checkpoint_info = (0.0, Checkpoints(), None, tmp_path)
+        mmp_client_object.get_checkpoint_info.return_value = checkpoint_info
         mmp_client.return_value = mmp_client_object
         ports = {
             Operator.F_INIT: ['in'],
@@ -159,9 +165,10 @@ def test_is_vector_port(instance):
 
 
 def test_send(instance, message):
+    instance._trigger_manager._cpts_considered_until = 17.0
     instance.send('out', message, 1)
     assert instance._communicator.send_message.called_with(
-            'out', message, 1)
+            'out', message, 1, 17.0)
 
 
 def test_send_invalid_port(instance, message):
@@ -232,7 +239,8 @@ def test_reuse_instance_receive_overlay(instance):
     test_overlay = Settings()
     test_overlay['test2'] = 'abc'
     recv = instance._communicator.receive_message
-    recv.return_value = Message(0.0, None, test_overlay, test_base_settings)
+    msg = Message(0.0, None, test_overlay, test_base_settings)
+    recv.return_value = msg, 0.0
     instance.reuse_instance()
     assert instance._communicator.receive_message.called_with(
         'muscle_settings_in')
@@ -244,9 +252,9 @@ def test_reuse_instance_receive_overlay(instance):
 def test_reuse_instance_closed_port(instance):
     def receive_message(port_name, slot=None, default=None):
         if port_name == 'muscle_settings_in':
-            return Message(0.0, None, Settings(), Settings())
+            return Message(0.0, None, Settings(), Settings()), 0.0
         elif port_name == 'in':
-            return Message(0.0, None, ClosePort(), Settings())
+            return Message(0.0, None, ClosePort(), Settings()), 1.0
         assert False    # pragma: no cover
 
     def get_port(port_name):
@@ -272,10 +280,10 @@ def get_port(port_name):
 def test_reuse_instance_vector_port(instance2):
     def receive_message(port_name, slot=None, default=None):
         if port_name == 'muscle_settings_in':
-            return Message(0.0, None, Settings(), Settings())
+            return Message(0.0, None, Settings(), Settings()), 0.0
         elif port_name == 'in':
             data = 'test {}'.format(slot)
-            return Message(0.0, None, data, Settings())
+            return Message(0.0, None, data, Settings()), 0.0
         assert False    # pragma: no cover
 
     instance2._communicator.receive_message = receive_message
@@ -300,7 +308,7 @@ def receive_message(port_name, slot=None, default=None):
 
 def test_reuse_instance_no_f_init_ports(instance):
     instance._communicator.receive_message.return_value = Message(
-            0.0, None, Settings(), Settings())
+            0.0, None, Settings(), Settings()), 0.0
     instance._communicator.list_ports.return_value = {}
     instance._communicator.settings_in_connected.return_value = False
     do_reuse = instance.reuse_instance()
diff --git a/libmuscle/python/libmuscle/test/test_mpp_message.py b/libmuscle/python/libmuscle/test/test_mpp_message.py
index 79ee8ee6..aaebe351 100644
--- a/libmuscle/python/libmuscle/test/test_mpp_message.py
+++ b/libmuscle/python/libmuscle/test/test_mpp_message.py
@@ -15,16 +15,20 @@ def test_create() -> None:
     timestamp = 10.0
     next_timestamp = 11.0
     settings_overlay = (6789).to_bytes(2, 'little', signed=True)
+    message_number = 0
+    saved_until = 1.6
     data = (12345).to_bytes(2, 'little', signed=True)
     msg = MPPMessage(
             sender, receiver, None, timestamp, next_timestamp,
-            settings_overlay, data)
+            settings_overlay, message_number, saved_until, data)
     assert msg.sender == sender
     assert msg.receiver == receiver
     assert msg.port_length is None
     assert msg.timestamp == 10.0
     assert msg.next_timestamp == 11.0
     assert msg.settings_overlay == settings_overlay
+    assert msg.message_number == 0
+    assert msg.saved_until == 1.6
     assert msg.data == data
 
 
@@ -43,7 +47,7 @@ def test_grid_encode() -> None:
     grid = Grid(array, ['x', 'y', 'z'])
     msg = MPPMessage(
             sender, receiver, None, timestamp, next_timestamp, Settings(),
-            grid)
+            0, 1.0, grid)
 
     wire_data = msg.encoded()
     mcp_decoded = msgpack.unpackb(wire_data, raw=False)
@@ -86,6 +90,8 @@ def test_grid_decode() -> None:
             'timestamp': 0.0,
             'next_timestamp': None,
             'settings_overlay': msgpack.ExtType(1, settings_data),
+            'message_number': 0,
+            'saved_until': 9.9,
             'data': msgpack.ExtType(2, grid_data)}
 
     wire_data = msgpack.packb(msg_dict, use_bin_type=True)
@@ -135,7 +141,7 @@ def test_grid_roundtrip() -> None:
         grid = Grid(array, ['x', 'y', 'z'])
         msg = MPPMessage(
                 sender, receiver, None, timestamp, next_timestamp, Settings(),
-                grid)
+                0, 1.0, grid)
 
         wire_data = msg.encoded()
         msg_out = MPPMessage.from_bytes(wire_data)
@@ -169,7 +175,7 @@ def test_non_contiguous_grid_roundtrip() -> None:
     grid = Grid(array.real, ['a', 'b', 'c'])
     msg = MPPMessage(
             sender, receiver, None, timestamp, next_timestamp, Settings(),
-            grid)
+            0, 7.7, grid)
 
     wire_data = msg.encoded()
     msg_out = MPPMessage.from_bytes(wire_data)
diff --git a/libmuscle/python/libmuscle/test/test_outbox.py b/libmuscle/python/libmuscle/test/test_outbox.py
index 6b22f068..a2e97c40 100644
--- a/libmuscle/python/libmuscle/test/test_outbox.py
+++ b/libmuscle/python/libmuscle/test/test_outbox.py
@@ -19,6 +19,7 @@ def message():
             Ref('sender.out'), Ref('receiver.in'),
             None, 0.0, 1.0,
             bytes(),
+            0, 1.0,
             'testing'.encode('utf-8'))
 
 
diff --git a/libmuscle/python/libmuscle/test/test_snapshot.py b/libmuscle/python/libmuscle/test/test_snapshot.py
new file mode 100644
index 00000000..f459a001
--- /dev/null
+++ b/libmuscle/python/libmuscle/test/test_snapshot.py
@@ -0,0 +1,75 @@
+import pytest
+from ymmsl import Settings
+
+from libmuscle.communicator import Message
+from libmuscle.snapshot import Snapshot, MsgPackSnapshot, SnapshotMetadata
+
+
+@pytest.fixture
+def snapshot() -> Snapshot:
+    triggers = ['test triggers']
+    wallclock_time = 15.3
+    port_message_counts = {'in': [1], 'out': [4], 'muscle_settings_in': [0]}
+    is_final = True
+    message = Message(1.2, None, 'test_data')
+    snapshot = MsgPackSnapshot(
+            triggers, wallclock_time, port_message_counts, is_final, message)
+    assert snapshot.triggers == triggers
+    assert snapshot.wallclock_time == wallclock_time
+    assert snapshot.port_message_counts == port_message_counts
+    assert snapshot.is_final_snapshot == is_final
+    assert snapshot.message == message
+    return snapshot
+
+
+def test_snapshot(snapshot: Snapshot) -> None:
+    assert isinstance(snapshot, Snapshot)
+
+    binary_snapshot = snapshot.to_bytes()
+    assert isinstance(binary_snapshot, bytes)
+
+    snapshot2 = MsgPackSnapshot.from_bytes(binary_snapshot)
+
+    assert snapshot2.triggers == snapshot.triggers
+    assert snapshot2.wallclock_time == snapshot.wallclock_time
+    assert snapshot2.port_message_counts == snapshot.port_message_counts
+    assert snapshot2.is_final_snapshot == snapshot.is_final_snapshot
+    assert snapshot2.message.timestamp == snapshot.message.timestamp
+    assert snapshot2.message.next_timestamp == snapshot.message.next_timestamp
+    assert snapshot2.message.data == snapshot.message.data
+
+
+def test_snapshot_metadata(snapshot: Snapshot) -> None:
+    metadata = SnapshotMetadata.from_snapshot(snapshot, 'test')
+
+    assert metadata.triggers == snapshot.triggers
+    assert metadata.wallclock_time == snapshot.wallclock_time
+    assert metadata.port_message_counts == snapshot.port_message_counts
+    assert metadata.is_final_snapshot == snapshot.is_final_snapshot
+    assert metadata.timestamp == snapshot.message.timestamp
+    assert metadata.next_timestamp == snapshot.message.next_timestamp
+    assert metadata.snapshot_filename == 'test'
+
+
+def test_message_with_settings() -> None:
+    message = Message(1.0, 2.0, 'test_data', Settings({'setting': True}))
+    snapshot = MsgPackSnapshot([], 0, {}, False, message)
+    assert snapshot.message.settings.get('setting') is True
+
+    binary_snapshot = snapshot.to_bytes()
+    assert isinstance(binary_snapshot, bytes)
+
+    snapshot2 = MsgPackSnapshot.from_bytes(binary_snapshot)
+    assert snapshot2.message.settings.get('setting') is True
+
+
+def test_implicit_snapshot() -> None:
+    message = None
+    snapshot = MsgPackSnapshot([], 0, {}, True, message)
+    assert snapshot.message is None
+
+    binary_snapshot = snapshot.to_bytes()
+    assert isinstance(binary_snapshot, bytes)
+
+    snapshot2 = MsgPackSnapshot.from_bytes(binary_snapshot)
+    assert snapshot2.message is None
diff --git a/libmuscle/python/libmuscle/test/test_snapshot_manager.py b/libmuscle/python/libmuscle/test/test_snapshot_manager.py
new file mode 100644
index 00000000..31423bb0
--- /dev/null
+++ b/libmuscle/python/libmuscle/test/test_snapshot_manager.py
@@ -0,0 +1,121 @@
+from pathlib import Path
+from unittest.mock import MagicMock
+
+from ymmsl import Reference
+
+from libmuscle.communicator import Message
+from libmuscle.snapshot import SnapshotMetadata
+from libmuscle.snapshot_manager import SnapshotManager
+
+
+def test_no_checkpointing(tmp_path: Path) -> None:
+    manager = MagicMock()
+    communicator = MagicMock()
+    communicator.get_message_counts.return_value = {}
+    snapshot_manager = SnapshotManager(Reference('test'), manager, communicator)
+
+    snapshot_manager.prepare_resume(None, tmp_path)
+    assert not snapshot_manager.resuming_from_intermediate()
+    assert not snapshot_manager.resuming_from_final()
+
+
+def test_save_load_snapshot(tmp_path: Path) -> None:
+    manager = MagicMock()
+    communicator = MagicMock()
+    port_message_counts = {'in': [1], 'out': [2], 'muscle_settings_in': [0]}
+    communicator.get_message_counts.return_value = port_message_counts
+
+    instance_id = Reference('test[1]')
+    snapshot_manager = SnapshotManager(instance_id, manager, communicator)
+
+    snapshot_manager.prepare_resume(None, tmp_path)
+    assert not snapshot_manager.resuming_from_intermediate()
+    assert not snapshot_manager.resuming_from_final()
+
+    snapshot_manager.save_snapshot(
+            Message(0.2, None, 'test data'), False, ['test'], 13.0)
+
+    communicator.get_message_counts.assert_called_with()
+    manager.submit_snapshot_metadata.assert_called()
+    instance, metadata = manager.submit_snapshot_metadata.call_args[0]
+    assert instance == instance_id
+    assert isinstance(metadata, SnapshotMetadata)
+    assert metadata.triggers == ['test']
+    assert metadata.wallclock_time == 13.0
+    assert metadata.timestamp == 0.2
+    assert metadata.next_timestamp is None
+    assert metadata.port_message_counts == port_message_counts
+    assert not metadata.is_final_snapshot
+    snapshot_path = Path(metadata.snapshot_filename)
+    assert snapshot_path.parent == tmp_path
+    assert snapshot_path.name == 'test-1_1.pack'
+
+    snapshot_manager2 = SnapshotManager(instance_id, manager, communicator)
+
+    snapshot_manager2.prepare_resume(snapshot_path, tmp_path)
+    communicator.restore_message_counts.assert_called_with(port_message_counts)
+
+    assert snapshot_manager2.resuming_from_intermediate()
+    assert not snapshot_manager2.resuming_from_final()
+    msg = snapshot_manager2.load_snapshot()
+    assert msg.timestamp == 0.2
+    assert msg.next_timestamp is None
+    assert msg.data == 'test data'
+
+    snapshot_manager2.save_snapshot(
+            Message(0.6, None, 'test data2'), True, ['test'], 42.2, 1.2)
+
+    instance, metadata = manager.submit_snapshot_metadata.call_args[0]
+    assert instance == instance_id
+    assert isinstance(metadata, SnapshotMetadata)
+    assert metadata.triggers == ['test']
+    assert metadata.wallclock_time == 42.2
+    assert metadata.timestamp == 0.6
+    assert metadata.next_timestamp is None
+    assert metadata.port_message_counts == port_message_counts
+    assert metadata.is_final_snapshot
+    snapshot_path = Path(metadata.snapshot_filename)
+    assert snapshot_path.parent == tmp_path
+    assert snapshot_path.name == 'test-1_3.pack'
+
+    assert snapshot_manager2.resuming_from_intermediate()
+    assert not snapshot_manager2.resuming_from_final()
+    snapshot_manager2.load_snapshot()
+    assert snapshot_manager2.resuming_from_intermediate()
+    assert not snapshot_manager2.resuming_from_final()
+
+
+def test_save_load_implicit_snapshot(tmp_path: Path) -> None:
+    manager = MagicMock()
+    communicator = MagicMock()
+    port_message_counts = {'in': [1], 'out': [2], 'muscle_settings_in': [0]}
+    communicator.get_message_counts.return_value = port_message_counts
+
+    instance_id = Reference('test[1]')
+    snapshot_manager = SnapshotManager(instance_id, manager, communicator)
+
+    snapshot_manager.prepare_resume(None, tmp_path)
+
+    assert not snapshot_manager.resuming_from_intermediate()
+    assert not snapshot_manager.resuming_from_final()
+    # save implicit snapshot
+    snapshot_manager.save_snapshot(None, True, ['implicit'], 1.0, 1.5)
+
+    manager.submit_snapshot_metadata.assert_called_once()
+    instance, metadata = manager.submit_snapshot_metadata.call_args[0]
+    assert instance == instance_id
+    assert isinstance(metadata, SnapshotMetadata)
+    snapshot_path = Path(metadata.snapshot_filename)
+    manager.submit_snapshot_metadata.reset_mock()
+
+    snapshot_manager2 = SnapshotManager(instance_id, manager, communicator)
+
+    snapshot_manager2.prepare_resume(snapshot_path, tmp_path)
+    communicator.restore_message_counts.assert_called_with(port_message_counts)
+    manager.submit_snapshot_metadata.assert_called_once()
+    manager.submit_snapshot_metadata.reset_mock()
+
+    assert not snapshot_manager2.resuming_from_intermediate()
+    assert not snapshot_manager2.resuming_from_final()
+    snapshot_manager2.save_snapshot(None, True, ['implicit'], 12.3, 2.5)
+    manager.submit_snapshot_metadata.assert_called_once()