diff --git a/.gitignore b/.gitignore
index cde952be38688dc3206700e1a73902bafc46defb..7b710fb915652dba267801bff76a20cb9c063475 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1 +1,3 @@
-TODO.md
\ No newline at end of file
+TODO.md
+docs/planning/very_first_results.zip
+*results*/
\ No newline at end of file
diff --git a/src/Persephone.jl b/src/Persephone.jl
index ad1fd4f6ae2cd5ca2526d44a79c6abc8066b69fd..a8becf93d5b91193498441e630f29bc0717e6a2d 100644
--- a/src/Persephone.jl
+++ b/src/Persephone.jl
@@ -42,6 +42,7 @@ include("core/landscape.jl")
 include("farm/farm.jl")
 include("crop/crops.jl")
 include("nature/nature.jl")
+include("nature/ecologicaldata.jl")
 include("nature/wolpertinger.jl")
 include("nature/wyvern.jl")
 include("core/simulation.jl")
diff --git a/src/core/input.jl b/src/core/input.jl
index 229e016350c4a38755a761d213b81a4301cbbc17..2ba80a21522c3e26f714174cbb19d63adea81813 100644
--- a/src/core/input.jl
+++ b/src/core/input.jl
@@ -31,11 +31,11 @@ let settings::Dict{String, Dict{String, Any}}
     end
 
     """
-        astoml(io)
+        printparameters(io)
 
     Print all settings in TOML format to the given output stream.
     """
-    global function astoml(stream::IO=stdout)
+    global function printparameters(stream::IO=stdout)
         TOML.print(stream, settings)
     end
 end
diff --git a/src/core/output.jl b/src/core/output.jl
index eb767143ae1447c41b3ea323c8b248a6f44425a1..87afd9704eba6d577636e8e22391e93540708f33 100644
--- a/src/core/output.jl
+++ b/src/core/output.jl
@@ -4,8 +4,6 @@
 ###
 
 const LOGFILE = "simulation.log"
-const POPFILE = "populations.csv"
-const INDFILE = "individuals.csv"
 
 ## Note: `setupdatadir()` was adapted from the GeMM model by Leidinger et al.
 ## (https://github.com/CCTB-Ecomods/gemm/blob/master/src/output.jl)
@@ -46,7 +44,7 @@ function setupdatadir()
         println(f, "# This file was generated automatically.")
         println(f, "# Simulation run on $(string(Dates.format(Dates.now(), "d u Y HH:MM:SS"))),")
         println(f, "# with git commit $(read(`git rev-parse HEAD`, String))#\n")
-        astoml(f)
+        printparameters(f)
     end
     # Copy the map files to the output folder
     lcmap = param("core.landcovermap")
@@ -55,109 +53,80 @@ function setupdatadir()
     !(isfile(ffmap)) && Base.error("The map file $(ffmap) doesn't exist.")
     cp(lcmap, joinpath(param("core.outdir"), basename(lcmap)), force = true)
     cp(ffmap, joinpath(param("core.outdir"), basename(ffmap)), force = true)
-    # Create the data output file(s)
-    initnaturedatafiles()
 end
 
+"""
+    DataOutput
 
-## NATURE MODEL OUTPUT FUNCTIONS
+A struct for organising model output. This is designed for text-based data output
+that is updated more or less regularly (e.g. population data in csv files).
+Submodels can register their own output functions using `newdataoutput()`.
 
-#XXX Should these be moved to a separate file (src/nature/natureoutput.jl)?
-# At the moment I've kept them here so that all output functions are in the same place.
-let nextpopout = nothing, nextindout = nothing
-    """
-        initnaturedatafiles()
+Struct fields:
+    - filename: the name of the file to be created in the user-specified output directory
+    - header: a string to be written to the start of the file as it is initialised
+    - outputfunction: a function that takes a model object and returns a string to write to file
+    - frequency: how often to call the output function (daily/monthly/yearly/end/never)
+"""
+struct DataOutput
+    filename::String
+    header::String
+    outputfunction::Function
+    frequency::String
+end
 
-    Initialise the files needed for the nature output data and set the
-    date for the next data output, depending on the model settings.
-    """
-    global function initnaturedatafiles()
-        startdate = param("core.startdate")
-        popoutfreq = param("nature.popoutfreq")
-        indoutfreq = param("nature.indoutfreq")
-        if popoutfreq != "never"
-            open(joinpath(param("core.outdir"), POPFILE), "w") do f
-                println(f, "Date;Species;Abundance")
-            end
-            (popoutfreq == "daily") && (nextpopout = startdate)
-            (popoutfreq == "monthly") && (nextpopout = startdate+Month(1))
-            (popoutfreq == "yearly") && (nextpopout = startdate+Year(1))
-            (popoutfreq == "end") && (nextpopout = param("core.enddate"))
-        else nextpopout = startdate - Day(1) end
-        if indoutfreq != "never"
-            open(joinpath(param("core.outdir"), INDFILE), "w") do f
-                println(f, "Date;ID;X;Y;Species;Sex;Age;Energy")
-            end
-            (indoutfreq == "daily") && (nextindout = startdate)
-            (indoutfreq == "monthly") && (nextindout = startdate+Month(1))
-            (indoutfreq == "yearly") && (nextindout = startdate+Year(1))
-            (indoutfreq == "end") && (nextindout = param("core.enddate"))
-        else nextindout = startdate - Day(1) end
-    end
+let outputregistry = Vector{DataOutput}(),
+    nextmonthlyoutput = today(),
+    nextyearlyoutput = today()
 
     """
-        recordnaturedata(model)
+        newdataoutput(filename, header, outputfunction, frequency)
 
-    Record the relevant data output from the nature model, depending on the frequency settings.
+    Create and register a new data output. This function must be called by all submodels
+    that want to have their output functions called regularly.
     """
-    global function recordnaturedata(model::AgentBasedModel)
-        popoutfreq = param("nature.popoutfreq")
-        indoutfreq = param("nature.indoutfreq")
-        if model.date == nextpopout
-            savepopulationdata(model)
-            (popoutfreq == "daily") && (nextpopout = model.date+Day(1))
-            (popoutfreq == "monthly") && (nextpopout = model.date+Month(1))
-            (popoutfreq == "yearly") && (nextpopout = model.date+Year(1))
+    global function newdataoutput(filename::String, header::String,
+                                  outputfunction::Function, frequency::String)
+        if !(frequency in ("daily", "monthly", "yearly", "end", "never"))
+            Base.error("Invalid frequency '$frequency' for $filename.")
         end
-        if model.date == nextindout
-            saveindividualdata(model)
-            (indoutfreq == "daily") && (nextindout = model.date+Day(1))
-            (indoutfreq == "monthly") && (nextindout = model.date+Month(1))
-            (indoutfreq == "yearly") && (nextindout = model.date+Year(1))
+        ndo = DataOutput(filename, header, outputfunction, frequency)
+        append!(outputregistry, [ndo])
+        if frequency != "never"
+            open(joinpath(param("core.outdir"), filename), "w") do f
+                println(f, header)
+            end
         end
     end
-end
-
-##XXX The current method of controlling the next output date for the different
-## output functions is a bit messy. If the farm model introduces similar parameters
-## to popoutfreq/indoutfreq, it would be worth revising this to a cleaner hook system.
-
-"""
-    savepopulationdata(model)
 
-Print a comma-separated set of lines to `populations.csv`, giving the current date
-and population size for each animal species. May be called never, daily, monthly,
-yearly, or at the end of a simulation, depending on the parameter `nature.popoutfreq`.
-"""
-function savepopulationdata(model::AgentBasedModel)
-    pops = Dict{String,Int}(s=>0 for s = param("nature.targetspecies"))
-    for a in allagents(model)
-        (typeof(a) != Animal) && continue
-        pops[a.species.name] += 1
-    end
-    open(joinpath(param("core.outdir"), POPFILE), "a") do f
-        for p in keys(pops)
-            println(f, join([model.date, p, pops[p]], ";"))
+    """
+        outputdata(model)
+    
+    Cycle through all registered data outputs and activate them according to their
+    configured frequency.
+    """
+    global function outputdata(model::AgentBasedModel)
+        #XXX all output functions are run on the first update (regardless of frequency)
+        # -> should they all be run on the last update, too?
+        if model.date == param("core.startdate")
+            nextmonthlyoutput = model.date
+            nextyearlyoutput = model.date
         end
-    end
-end
-
-"""
-    saveindividualdata(model)
-
-Print a comma-separated set of lines to `individuals.csv`, listing all properties
-of all animal individuals in the model. May be called never, daily, monthly, yearly, or
-at the end of a simulation, depending on the parameter `nature.indoutfreq`.
-WARNING: Produces very big files!
-"""
-function saveindividualdata(model::AgentBasedModel)
-    data = ""
-    for a in allagents(model)
-        (typeof(a) != Animal) && continue
-        entry = join([model.date,a.id,a.pos[1],a.pos[2],a.species.name,a.sex,a.age,a.energy], ";")
-        data = data*entry*"\n"
-    end
-    open(joinpath(param("core.outdir"), INDFILE), "a") do f
-        print(f, data)
+        for output in outputregistry
+            (output.frequency == "never") && continue
+            # check if this output should be activated today
+            if (output.frequency == "daily") ||
+                (output.frequency == "monthly" && model.date == nextmonthlyoutput) ||
+                (output.frequency == "yearly" && model.date == nextyearlyoutput) ||
+                (output.frequency == "end" && model.date == param("core.enddate"))
+                open(joinpath(param("core.outdir"), output.filename), "a") do f
+                    outstring = output.outputfunction(model)
+                    (outstring[end] != '\n') && (outstring *= '\n')
+                    print(f, outstring)
+                end                
+            end
+        end
+        (model.date == nextmonthlyoutput) && (nextmonthlyoutput = model.date + Month(1))
+        (model.date == nextyearlyoutput) && (nextyearlyoutput = model.date + Year(1))
     end
 end
diff --git a/src/core/simulation.jl b/src/core/simulation.jl
index 36511a4e1e3f21ed0f6589ef75ef2a89b7466a6d..7ef064e39bbd24f3ad0c4b5b47d036a51b033a53 100644
--- a/src/core/simulation.jl
+++ b/src/core/simulation.jl
@@ -40,7 +40,7 @@ function stepsimulation!(model::AgentBasedModel)
         #The animal may have been killed, so we need a try/catch
         try stepagent!(getindex(model, a), model) catch keyerror end
     end
-    recordnaturedata(model)
+    outputdata(model)
     model.date += Day(1)
 end
 
diff --git a/src/nature/ecologicaldata.jl b/src/nature/ecologicaldata.jl
new file mode 100644
index 0000000000000000000000000000000000000000..f01ead1be6cffc293493c4012a61035369e07699
--- /dev/null
+++ b/src/nature/ecologicaldata.jl
@@ -0,0 +1,46 @@
+### Persephone - a socio-economic-ecological model of European agricultural landscapes.
+###
+### This file includes the functions for collecting and saving ecological output data.
+###
+
+const POPFILE = "populations.csv"
+const INDFILE = "individuals.csv"
+
+"""
+    savepopulationdata(model)
+
+Print a comma-separated set of lines to `populations.csv`, giving the current date
+and population size for each animal species. May be called never, daily, monthly,
+yearly, or at the end of a simulation, depending on the parameter `nature.popoutfreq`.
+"""
+function savepopulationdata(model::AgentBasedModel)
+    pops = Dict{String,Int}(s=>0 for s = param("nature.targetspecies"))
+    for a in allagents(model)
+        (typeof(a) != Animal) && continue
+        pops[a.species.name] += 1
+    end
+    data = ""
+    for p in keys(pops)
+        data *= join([model.date, p, pops[p]], ";")*"\n"
+    end
+    data
+end
+
+"""
+    saveindividualdata(model)
+
+Print a comma-separated set of lines to `individuals.csv`, listing all properties
+of all animal individuals in the model. May be called never, daily, monthly, yearly, or
+at the end of a simulation, depending on the parameter `nature.indoutfreq`.
+WARNING: Produces very big files!
+"""
+function saveindividualdata(model::AgentBasedModel)
+    data = ""
+    for a in allagents(model)
+        (typeof(a) != Animal) && continue
+        entry = join([model.date,a.id,a.pos[1],a.pos[2],a.species.name,a.sex,a.age,a.energy], ";")
+        data *= entry*"\n"
+    end
+    data
+end
+
diff --git a/src/nature/nature.jl b/src/nature/nature.jl
index c41af8fc2d111b8a3a77545a584b63eca8195a98..3d6db6f2beb347436624ec460f6dd0636173e072 100644
--- a/src/nature/nature.jl
+++ b/src/nature/nature.jl
@@ -97,8 +97,13 @@ end
 Initialise the model with all simulated animal populations.
 """
 function initnature!(model::AgentBasedModel)
-    #The config file determines which species are simulated in this run
+    # The config file determines which species are simulated in this run
     for s in param("nature.targetspecies")
         getspecies(s).initpop!(model)
     end
+    # Initialise the data output
+    newdataoutput(POPFILE, "Date;Species;Abundance", savepopulationdata,
+                  param("nature.popoutfreq"))
+    newdataoutput(INDFILE, "Date;ID;X;Y;Species;Sex;Age;Energy",
+                  saveindividualdata, param("nature.indoutfreq"))
 end
diff --git a/test.sh b/test.sh
index 15fb6c9c136143ad1c8340a769504892c4b42498..8641b43fb7a920882e2b8180427b801cd272e81f 100755
--- a/test.sh
+++ b/test.sh
@@ -6,6 +6,6 @@ then
     rm -r test_results
 fi
 
-./run.jl -o test_results -s 1
+./run.jl -o test_results -s 2
 
 src/analysis/analyse_nature.R