diff --git a/src/PlantSimEngine.jl b/src/PlantSimEngine.jl
index a6ec6233f..ebf59dcb9 100644
--- a/src/PlantSimEngine.jl
+++ b/src/PlantSimEngine.jl
@@ -43,6 +43,9 @@ include("component_models/RefVector.jl")
 # Simulation table (time-step table, from PlantMeteo):
 include("component_models/TimeStepTable.jl")
 
+# Declaring the dependency graph
+include("dependencies/dependency_graph.jl")
+
 # List of models:
 include("component_models/ModelList.jl")
 include("mtg/MultiScaleModel.jl")
@@ -53,8 +56,7 @@ include("component_models/get_status.jl")
 # Transform into a dataframe:
 include("dataframe.jl")
 
-# Model dependencies:
-include("dependencies/dependency_graph.jl")
+# Computing model dependencies:
 include("dependencies/soft_dependencies.jl")
 include("dependencies/hard_dependencies.jl")
 include("dependencies/traversal.jl")
diff --git a/src/component_models/ModelList.jl b/src/component_models/ModelList.jl
index b0e23be6a..771a50e27 100644
--- a/src/component_models/ModelList.jl
+++ b/src/component_models/ModelList.jl
@@ -147,7 +147,8 @@ julia> [typeof(models[i][1]) for i in keys(status(models))]
 struct ModelList{M<:NamedTuple,S}
     models::M
     status::S
-    type_promotion::Union{Nothing, Dict}
+    type_promotion::Union{Nothing,Dict}
+    dependency_graph::DependencyGraph
 end
 
 #=function ModelList(models::M, status::Status) where {M<:NamedTuple{names,T} where {names,T<:NTuple{N,<:AbstractModel} where {N}}}
@@ -160,7 +161,6 @@ function ModelList(
     status=nothing,
     type_promotion::Union{Nothing,Dict}=nothing,
     variables_check::Bool=true,
-    nsteps=nothing,
     kwargs...
 )
 
@@ -187,10 +187,13 @@ function ModelList(
     ts_kwargs = homogeneous_ts_kwargs(status)
     ts_kwargs = add_model_vars(ts_kwargs, mods, type_promotion)
 
+
+
     model_list = ModelList(
         mods,
         ts_kwargs,
-        type_promotion
+        type_promotion,
+        dep(; verbose=true, mods...)
     )
     variables_check && !is_initialized(model_list)
 
@@ -219,8 +222,8 @@ function add_model_vars(x, models, type_promotion)
 
     # If the user gave a status, we check if all the variables are already initialized:
     vars_in_x = status_keys(x)
-    status_x = 
-    all([k in vars_in_x for k in keys(ref_vars)]) && return isa(x, Status) ? x : Status(x)  # If so, we return the input
+    status_x =
+        all([k in vars_in_x for k in keys(ref_vars)]) && return isa(x, Status) ? x : Status(x)  # If so, we return the input
 
     # Else, we add the variables by making a new object (carefull, this is a copy so it takes more time):
 
@@ -229,15 +232,15 @@ function add_model_vars(x, models, type_promotion)
 
     # If the user gave an empty status, we initialize all variables to their default values:
     if x === nothing
-       return Status(ref_vars)
+        return Status(ref_vars)
     end
-    
+
     if Tables.istable(x)
         # This situation only occurs if the user provided a table instead of a status
         # Meaning we have a status of vector values, all initialized up to a certain point
         # Unsure this is desirable, as that means run! does nothing or overwrites everything
         # Anyway, we wish to create a NamedTuple() of Vectors here
-        x_full = (;zip(propertynames(x), Tables.columns(x))...)
+        x_full = (; zip(propertynames(x), Tables.columns(x))...)
         x_full = merge(ref_vars, x_full)
 
     else
@@ -286,7 +289,7 @@ PlantSimEngine.homogeneous_ts_kwargs((Tₗ=[25.0, 26.0], aPPFD=1000.0))
 function homogeneous_ts_kwargs(kwargs::NamedTuple{N,T}) where {N,T}
     length(kwargs) == 0 && return kwargs
     vars_vals = collect(Any, values(kwargs))
- 
+
     vars_array = NamedTuple{keys(kwargs)}(j for j in vars_vals)
 
     return vars_array
@@ -325,7 +328,8 @@ function Base.copy(m::T) where {T<:ModelList}
     ModelList(
         m.models,
         deepcopy(m.status),
-        deepcopy(m.type_promotion)
+        deepcopy(m.type_promotion),
+        deepcopy(m.dependency_graph)
     )
 end
 
@@ -333,7 +337,8 @@ function Base.copy(m::T, status) where {T<:ModelList}
     ModelList(
         m.models,
         status,
-        deepcopy(m.type_promotion)
+        deepcopy(m.type_promotion),
+        deepcopy(m.dependency_graph)
     )
 end
 
@@ -465,7 +470,7 @@ function convert_vars!(mapped_vars::Dict{String,Dict{Symbol,Any}}, type_promotio
 end
 
 function Base.show(io::IO, ::MIME"text/plain", t::ModelList)
-    print(io, dep(t, verbose=false))
+    print(io, dep(t))
     print(io, status(t))
 end
 
diff --git a/src/component_models/Status.jl b/src/component_models/Status.jl
index 3294f3222..60cc07ad5 100644
--- a/src/component_models/Status.jl
+++ b/src/component_models/Status.jl
@@ -135,26 +135,28 @@ end
 
 # Returns a status with all vector variables replaced with their first value (ie a Status ready for simulation)
 # also returns a tuple of symbols corresponding to the vector variables
-function flatten_status(s::Status)
-    status_values_flattened = NamedTuple()
-    vector_variables = NamedTuple()
-
-    for (var, value) in zip(keys(s), s)
-        if length(value) > 1
-            vector_variables = (vector_variables..., var)
-            status_values_flattened = (status_values_flattened..., value[1])
-        else
-            status_values_flattened = (status_values_flattened..., value)
-        end
+function flatten_status(s::Status{T}) where {T}
+    n_vars_several_values = findall(x -> length(x) > 1, s)
+    if length(n_vars_several_values) == 0
+        return s, n_vars_several_values
+    else
+        return Status{keys(s)}(first.(values(s))), n_vars_several_values
     end
-
-    return Status(; zip(keys(s), status_values_flattened)...), vector_variables
 end
 
-# Update to the next timestep the variables that were passed in as vectors by the user
-function update_vector_variables(s::Status, sf::Status, vector_variables, i)
-    for vec in vector_variables
-        sf[vec] = s[vec][i]
+"""
+    set_variables_at_timestep!(status_timestep::Status, user_status::Status, variables_to_update, timestep)
+
+Update `status_timestep` to the current values at the `timestep` for all `variables_to_update` in the status provided by the user (`user_status`).
+The variables to update are given in `variables_to_update`, which is a vector of symbols.
+
+`status_timestep` is a status representing a single time-step. `user_status` is the status provided that gives values for variables that are not computed by any model.
+It may give constant values or vectors of values, in which case the `timestep` is used to select the value to use for the current time step.
+
+"""
+function set_variables_at_timestep!(status_timestep::Status, user_status::Status, variables_to_update, timestep)
+    for vec in variables_to_update
+        status_timestep[vec] = user_status[vec][timestep]
     end
 end
 
diff --git a/src/dependencies/dependencies.jl b/src/dependencies/dependencies.jl
index ca50e010e..9252e7b07 100644
--- a/src/dependencies/dependencies.jl
+++ b/src/dependencies/dependencies.jl
@@ -1,8 +1,9 @@
 dep(::T, nsteps=1) where {T<:AbstractModel} = NamedTuple()
 
 """
-    dep(m::ModelList, nsteps=1; verbose::Bool=true)
+    dep(m::ModelList)
     dep(mapping::Dict{String,T}; verbose=true)
+    dep!(m::ModelList, nsteps=1)
 
 Get the model dependency graph given a ModelList or a multiscale model mapping. If one graph is returned, 
 then all models are coupled. If several graphs are returned, then only the models inside each graph are coupled, and
@@ -34,7 +35,12 @@ to other scales if needed. Then we transform all these nodes into soft dependenc
 Then we traverse all these and we set nodes that need outputs from other nodes as inputs as children/parents.
 If a node has no dependency, it is set as a root node and pushed into a new Dict (independant_process_root). This Dict is the returned dependency graph. And 
 it presents root nodes as independent starting points for the sub-graphs, which are the models that are coupled together. We can then traverse each of 
-these graphs independently to retrieve the models that are coupled together, in the right order of execution.
+these graphs independently to r
+
+# Notes
+
+The difference between `dep(m::ModelList)` and `dep!(m::ModelList, nsteps)` is that the first one returns the dependency graph found in the model list, while the 
+second one returns the dependency graph with the specified number of steps, modifying the simulation IDs of each node in the graph (`simulation_id=fill(0, nsteps)`).
 
 # Examples
 
@@ -75,8 +81,18 @@ function dep(nsteps=1; verbose::Bool=true, vars...)
     return deps
 end
 
-function dep(m::ModelList, nsteps=1; verbose::Bool=true)
-    dep(nsteps; verbose=verbose, m.models...)
+function dep(m::ModelList)
+    m.dependency_graph
+end
+
+function dep!(m::ModelList, nsteps=1)
+    traverse_dependency_graph!(m.dependency_graph; visit_hard_dep=false) do node
+        if length(node.simulation_id) != nsteps
+            node.simulation_id = fill(0, nsteps)
+        end
+    end
+
+    return m.dependency_graph
 end
 
 
diff --git a/src/mtg/save_results.jl b/src/mtg/save_results.jl
index 67b569148..51baffcbd 100644
--- a/src/mtg/save_results.jl
+++ b/src/mtg/save_results.jl
@@ -112,7 +112,7 @@ julia> collect(keys(preallocated_vars["Leaf"]))
 
 function pre_allocate_outputs(statuses, statuses_template, reverse_multiscale_mapping, vars_need_init, outs, nsteps; type_promotion=nothing, check=true)
     outs_ = Dict{String,Vector{Symbol}}()
-    
+
     # default behaviour : track everything
     if isnothing(outs)
         for organ in keys(statuses)
@@ -130,7 +130,7 @@ function pre_allocate_outputs(statuses, statuses_template, reverse_multiscale_ma
         end
     end
 
-    len = Dict{String, Int}()
+    len = Dict{String,Int}()
     for (organ, vals) in outs_
         len[organ] = length(outs_[organ])
         unique!(outs_[organ])
@@ -210,7 +210,7 @@ function pre_allocate_outputs(statuses, statuses_template, reverse_multiscale_ma
     node_type = only(node_type)
 
     # I don't know if this function barrier is necessary
-    preallocated_outputs = Dict{String, Vector}()
+    preallocated_outputs = Dict{String,Vector}()
     complete_preallocation_from_types!(preallocated_outputs, nsteps, outs_, node_type, statuses_template)
     return preallocated_outputs
 end
@@ -218,7 +218,7 @@ end
 function complete_preallocation_from_types!(preallocated_outputs, nsteps, outs_, node_type, statuses_template)
     types = Vector{DataType}()
     for organ in keys(outs_)
-        
+
         outs_no_node = filter(x -> x != :node, outs_[organ])
 
         #types = [typeof(status_from_template(statuses_template[organ])[var]) for var in outs[organ]]
@@ -230,14 +230,14 @@ function complete_preallocation_from_types!(preallocated_outputs, nsteps, outs_,
         symbols_tuple = (:timestep, :node, outs_no_node...,)
         # using node_type.parameters[1] is clunky, but covers both NodeMTG and AbstractNodeMTG types
         values_tuple = (1, MultiScaleTreeGraph.Node((node_type.parameters[1])("/", "Uninitialized", 0, 0),), values...,)
-    
+
         # Dummy value to make accessing the type easier 
         # (empty arrays don't have references to an instance, so their types can't be inspected and manipulated as easily)
-        dummy_status = (;zip(symbols_tuple, values_tuple)...)
+        dummy_status = (; zip(symbols_tuple, values_tuple)...)
         data = typeof(Status(dummy_status))[]
         resize!(data, nsteps)
-        
-        for ii in 1:nsteps 
+
+        for ii in 1:nsteps
             data[ii] = Status(dummy_status)
         end
         preallocated_outputs[organ] = data
@@ -278,22 +278,22 @@ function save_results!(object::GraphSimulation, i)
         # So there may be possible simplifications (maybe no need for a function barrier, perhaps the resizing could be made a one-liner...)
         # But this should work without causing visible performance regressions on XPalm
         len = length(outs[organ])
-        if length(statuses[organ]) + index - 1 > len            
+        if length(statuses[organ]) + index - 1 > len
             min_required = max(length(statuses[organ]) + index - len, index)
-            
-            extra_length = 2*min_required - len
+
+            extra_length = 2 * min_required - len
             data = eltype(outs[organ])[]
             resize!(data, extra_length)
             dummy_value = NamedTuple(outs[organ][1])
             # TODO set timestep to 0 for clarity ?
-            
+
             # Using fill! caused Ref issues, so call a Status constructor here instead of passing a prebuilt value
             # This will avoid having all array entries point to the same ref but keep construction cost at a minimum
             for new_entry in 1:extra_length
                 data[new_entry] = Status(dummy_value)
             end
 
-            outs[organ] = cat(outs[organ], data, dims=1)           
+            outs[organ] = cat(outs[organ], data, dims=1)
             #println("len : ", len, " statuses #", length(statuses[organ]), " index ", index)
             #println("min_required : ", min_required, " extra_length ", extra_length, " new len ", length(outs[organ]))
         end
@@ -316,15 +316,9 @@ function copy_tracked_outputs_into_vector!(outs_organ, i, statuses_organ, tracke
     return j
 end
 
-
 function pre_allocate_outputs(m::ModelList, outs, nsteps; type_promotion=nothing, check=true)
-    
-    # NOTE : init_variables recreates a DependencyGraph, it's not great
-    # TODO : copy ?
-    out_vars_pre_type_promotion = merge(init_variables(m; verbose=false)...)
-    
-    # bit hacky, could be cleaned up
-    out_vars_all = convert_vars(out_vars_pre_type_promotion, m.type_promotion)
+    st, = flatten_status(status(m))
+    out_vars_all = convert_vars(st, type_promotion)
 
     out_keys_requested = Symbol[]
     if !isnothing(outs)
@@ -337,9 +331,10 @@ function pre_allocate_outputs(m::ModelList, outs, nsteps; type_promotion=nothing
 
     # default implicit behaviour, track everything
     if isempty(out_keys_requested)
-        out_vars_requested = out_vars_all
+        # We already have the status here, just repeating its value:
+        out_vars_requested = NamedTuple(out_vars_all)
     else
-        unexpected_outputs = setdiff(out_keys_requested, status_keys(status(m)))
+        unexpected_outputs = setdiff(out_keys_requested, keys(st))
 
         if !isempty(unexpected_outputs)
             e = string(
@@ -354,22 +349,21 @@ function pre_allocate_outputs(m::ModelList, outs, nsteps; type_promotion=nothing
                 @info e
                 [delete!(unexpected_outputs, i) for i in unexpected_outputs]
             end
-        end    
+        end
 
         out_defaults_requested = (out_vars_all[i] for i in out_keys_requested)
-        out_vars_requested = (;zip(out_keys_requested, out_defaults_requested)...)
+        out_vars_requested = (; zip(out_keys_requested, out_defaults_requested)...)
     end
 
-    outputs_timestep = fill(out_vars_requested, nsteps)
-    return TimeStepTable([Status(i) for i in outputs_timestep])
+    return TimeStepTable([Status(out_vars_requested) for i in Base.OneTo(nsteps)])
 end
 
 function save_results!(status_flattened::Status, outputs, i)
-    if length(outputs) == 0 
-        return 
+    if length(outputs) == 0
+        return
     end
     outs = outputs[i]
-    
+
     for var in keys(outs)
         outs[var] = status_flattened[var]
     end
diff --git a/src/processes/model_initialisation.jl b/src/processes/model_initialisation.jl
index 5da66bdc5..95a8cc762 100755
--- a/src/processes/model_initialisation.jl
+++ b/src/processes/model_initialisation.jl
@@ -65,8 +65,8 @@ mapping = Dict(
 to_initialize(mapping)
 ```
 """
-function to_initialize(m::ModelList; verbose::Bool=true)
-    needed_variables = to_initialize(dep(m; verbose=verbose))
+function to_initialize(m::ModelList)
+    needed_variables = to_initialize(dep(m))
     to_init = Dict{Symbol,Tuple}()
     for (process, vars) in needed_variables
         # default_values = needed_variables[:process1]
@@ -245,7 +245,7 @@ function init_variables(model::T; verbose::Bool=true) where {T<:AbstractModel}
 end
 
 function init_variables(m::ModelList; verbose::Bool=true)
-    init_variables(dep(m; verbose=verbose))
+    init_variables(dep(m))
 end
 
 function init_variables(m::DependencyGraph)
@@ -300,7 +300,7 @@ is_initialized(models)
 ```
 """
 function is_initialized(m::T; verbose=true) where {T<:ModelList}
-    var_names = to_initialize(m; verbose=verbose)
+    var_names = to_initialize(m)
 
     if any([length(to_init) > 0 for (process, to_init) in pairs(var_names)])
         verbose && @info "Some variables must be initialized before simulation: $var_names (see `to_initialize()`)" maxlog = 1
diff --git a/src/run.jl b/src/run.jl
index 8958127a3..c0fd916a5 100644
--- a/src/run.jl
+++ b/src/run.jl
@@ -91,7 +91,7 @@ run!
 function adjust_weather_timesteps_to_given_length(desired_length, meteo)
     # This isn't ideal in terms of codeflow, but check_dimensions will kick in later
     # And determine whether there is a status vector length discrepancy
-     
+
     meteo_adjusted = meteo
 
     if DataFormat(meteo_adjusted) == TableAlike()
@@ -151,7 +151,7 @@ function run!(
 ) where {T<:Union{AbstractArray,AbstractDict},A}
 
     if executor != SequentialEx()
-         @warn string(
+        @warn string(
             "Parallelisation over objects was removed, (but may be reintroduced in the future). Parallelisation will only occur over timesteps."
         ) maxlog = 1
     end
@@ -161,7 +161,7 @@ function run!(
     # Each object:
     for obj in object
 
-        if isa(object, AbstractArray) 
+        if isa(object, AbstractArray)
             push!(outputs_collection, run!(obj, meteo, constants, extra, tracked_outputs=tracked_outputs, check=check, executor=executor))
         else
             outputs_collection[obj.first] = run!(obj.second, meteo, constants, extra, tracked_outputs=tracked_outputs, check=check, executor=executor)
@@ -183,11 +183,11 @@ function run!(
     check=true,
     executor=ThreadedEx()
 ) where {T<:ModelList}
-    
+
     meteo_adjusted = adjust_weather_timesteps_to_given_length(get_status_vector_max_length(object.status), meteo)
     nsteps = get_nsteps(meteo_adjusted)
 
-    dep_graph = dep(object, nsteps)
+    dep_graph = dep!(object, nsteps)
 
     if check
         # Check if the meteo data and the status have the same length (or length 1)
@@ -201,7 +201,7 @@ function run!(
         end
     end
 
-    
+
     if executor != SequentialEx() && nsteps > 1
         if !timestep_parallelizable(dep_graph)
             is_ts_parallel = which_timestep_parallelizable(dep_graph)
@@ -211,19 +211,19 @@ function run!(
                 "A parallel executor was provided (`executor=$(executor)`) but some models cannot be run in parallel: $mods_not_parallel. ",
                 "The simulation will be run sequentially. Use `executor=SequentialEx()` to remove this warning."
             ) maxlog = 1
-        else           
-            outputs_preallocated_mt = pre_allocate_outputs(object, tracked_outputs, nsteps)
+        else
+            outputs_preallocated_mt = pre_allocate_outputs(object, tracked_outputs, nsteps; type_promotion=object.type_promotion, check=check)
             local vars = length(outputs_preallocated_mt) > 0 ? keys(outputs_preallocated_mt[1]) : NamedTuple()
             status_flattened_template, vector_variables_mt = flatten_status(object.status)
 
             # Computing time-steps in parallel:
-            @floop executor for i in 1:nsteps                
-                @init begin 
+            @floop executor for i in 1:nsteps
+                @init begin
                     status_flattened = deepcopy(status_flattened_template)
                     roots = collect(dep_graph.roots)
                 end
                 meteo_i = meteo_adjusted[i]
-                update_vector_variables(object.status, status_flattened, vector_variables_mt, i)
+                set_variables_at_timestep!(status_flattened, status(object), vector_variables_mt, i)
                 for (process, node) in roots
                     run_node!(object, node, i, status_flattened, meteo_i, constants, extra)
                 end
@@ -235,8 +235,8 @@ function run!(
         end
     end
 
-    outputs_preallocated = pre_allocate_outputs(object, tracked_outputs, nsteps)
-    status_flattened, vector_variables = flatten_status(object.status)
+    outputs_preallocated = pre_allocate_outputs(object, tracked_outputs, nsteps; type_promotion=object.type_promotion, check=check)
+    status_flattened, vector_variables = flatten_status(status(object))
 
     # Not parallelizable over time-steps, it means some values depend on the previous value.
     # In this case we propagate the values of the variables from one time-step to the other, except for 
@@ -247,7 +247,7 @@ function run!(
     if nsteps == 1
         for (process, node) in roots
             run_node!(object, node, 1, status_flattened, meteo_adjusted, constants, extra)
-        end        
+        end
         save_results!(status_flattened, outputs_preallocated, 1)
     else
 
@@ -256,10 +256,10 @@ function run!(
                 run_node!(object, node, i, status_flattened, meteo_i, constants, extra)
             end
             save_results!(status_flattened, outputs_preallocated, i)
-            i + 1 <= nsteps && update_vector_variables(object.status, status_flattened, vector_variables, i + 1)
+            i + 1 <= nsteps && set_variables_at_timestep!(status_flattened, status(object), vector_variables, i + 1)
         end
     end
-    
+
     return outputs_preallocated
 end
 
@@ -273,7 +273,7 @@ function run!(
     tracked_outputs=nothing,
     check=true,
     executor=ThreadedEx()
-) where {T<:Union{AbstractArray, AbstractDict}}
+) where {T<:Union{AbstractArray,AbstractDict}}
 
     dep_graphs = [dep(obj) for obj in collect(values(object))]
     #obj_parallelizable = all([object_parallelizable(graph) for graph in dep_graphs])
@@ -305,7 +305,7 @@ function run!(
 
     # Each object:
     for obj in object
-        if isa(object, AbstractArray) 
+        if isa(object, AbstractArray)
             push!(outputs_collection, run!(obj, meteo, constants, extra, tracked_outputs=tracked_outputs, check=check, executor=executor))
         else
             outputs_collection[obj.first] = run!(obj.second, meteo, constants, extra, tracked_outputs=tracked_outputs, check=check, executor=executor)
diff --git a/test/test-ModelList.jl b/test/test-ModelList.jl
index b21df0bef..be62807b6 100644
--- a/test/test-ModelList.jl
+++ b/test/test-ModelList.jl
@@ -18,7 +18,6 @@
     leaf = ModelList(
         process1=Process1Model(1.0),
         process2=Process2Model(),
-        nsteps=3
     )
 
     @test length(status(leaf)) == 5
@@ -81,7 +80,6 @@ end;
         process1=Process1Model(1.0),
         process2=Process2Model(),
         status=(var1=15.0,),
-        nsteps=3
     )
 
     @test length(status(leaf)) == 5
@@ -123,8 +121,7 @@ end;
     inits = init_variables(leaf)
     sorted_vars = sort([keys(inits.process3)...])
 
-    @test [getfield(inits.process3, i) for i in sorted_vars] ==
-          fill(-Inf, 3)
+    @test [getfield(inits.process3, i) for i in sorted_vars] == fill(-Inf, 3)
 end;
 
 @testset "Copy a ModelList" begin
@@ -231,40 +228,38 @@ end=#
 @testset "ModelList outputs preallocation" begin
     meteo_day = CSV.read(joinpath(pkgdir(PlantSimEngine), "examples/meteo_day.csv"), DataFrame, header=18)
     vals = (var1=15.0, var2=0.3, TT_cu=cumsum(meteo_day.TT))
-    leaf =  ModelList(
+    leaf = ModelList(
         process1=Process1Model(1.0),
         process2=Process2Model(),
         status=vals
     )
-    outs=(:var3,)
+    outs = (:var3,)
 
     mtg, mapping, outputs_mapping, nsteps, filtered_outputs_modellist = test_filtered_output_begin(leaf, vals, outs, meteo_day)
     @test test_filtered_output(mtg, mapping, nsteps, outputs_mapping, meteo_day, filtered_outputs_modellist)
 
-    meteos = 
-    [Atmosphere(T=20.0, Wind=1.0, P=101.3, Rh=0.65, Ri_PAR_f=300.0),
-    CSV.read(joinpath(pkgdir(PlantSimEngine), "examples/meteo_day.csv"), DataFrame, header=18),
-    ]
+    meteos =
+        [Atmosphere(T=20.0, Wind=1.0, P=101.3, Rh=0.65, Ri_PAR_f=300.0),
+            CSV.read(joinpath(pkgdir(PlantSimEngine), "examples/meteo_day.csv"), DataFrame, header=18),
+        ]
     modellists, status_tuples, outs_vectors = get_modellist_bank()
 
     # remove some of the currently unhandled cases
-    outs_vectors = 
-    [
-        # this one has one tuple with a duplicate, and one with a nonexistent variable
-        [(:var1,), #=(:var1, :var1),=# (:var1, :var2), (:var1, :var3), (:var1, :var4, :var5), 
-        #=(:var2, :var7, :var3, :var1),=# (:var1, :var2, :var3, :var4, :var5)], 
-        [#=NamedTuple(),=# (:TT_cu,), (:TT_cu,:LAI) , (:biomass,:LAI), (:TT_cu, :LAI, :aPPFD, :biomass, :biomass_increment),], 
-        [#=NamedTuple(),=# (:var1,), (:var1, :var4), (:var1, :var2), (:var1, :var3), (:var1, :var4, :var6, :var5), 
-        #=(:var2, :var7, :var3, :var1),=# (:var1, :var2, :var3, :var4, :var5, :var6)], 
-        [#=NamedTuple(),=# (:var1,), (:var1, :var4), (:var1, :var2), (:var1, :var3), (:var1, :var4, :var6, :var5), 
-        (:var2, :var7, :var3, :var1), (:var1, :var2, :var3, :var4, :var5, :var6)],     
-        [#=NamedTuple(),=# (:var1,), (:var1, :var4), (:var1, :var2), (:var1, :var3), (:var1, :var4, :var6, :var5), 
-        (:var2, :var7, :var3, :var1), (:var1, :var2, :var3, :var4, :var5, :var6)
-        , (:var1, :var2, :var3, :var4, :var5, :var6, :var7, :var8, :var9)], 
-        [#=NamedTuple(),=# (:var1,), #=(:var1, :var1),=# (:var1, :var2), (:var1, :var3), (:var1, :var4, :var6, :var5), 
-        (:var2, :var7, :var3, :var1), (:var1, :var2, :var3, :var4, :var5, :var6)
-        , (:var1, :var2, :var3, :var4, :var5, :var6, :var7, #=:var8, :var9,=# :var0)], 
-    ]
+    outs_vectors =
+        [
+            # this one has one tuple with a duplicate, and one with a nonexistent variable
+            [(:var1,), (:var1, :var2), (:var1, :var3), (:var1, :var4, :var5), #=(:var1, :var1),=#
+                (:var1, :var2, :var3, :var4, :var5)],        #=(:var2, :var7, :var3, :var1),=#
+            [(:TT_cu,), (:TT_cu, :LAI), (:biomass, :LAI), (:TT_cu, :LAI, :aPPFD, :biomass, :biomass_increment),], #=NamedTuple(),=#
+            [(:var1,), (:var1, :var4), (:var1, :var2), (:var1, :var3), (:var1, :var4, :var6, :var5), #=NamedTuple(),=#
+                (:var1, :var2, :var3, :var4, :var5, :var6)],        #=(:var2, :var7, :var3, :var1),=#
+            [(:var1,), (:var1, :var4), (:var1, :var2), (:var1, :var3), (:var1, :var4, :var6, :var5), #=NamedTuple(),=#
+                (:var2, :var7, :var3, :var1), (:var1, :var2, :var3, :var4, :var5, :var6)],
+            [(:var1,), (:var1, :var4), (:var1, :var2), (:var1, :var3), (:var1, :var4, :var6, :var5), #=NamedTuple(),=#
+                (:var2, :var7, :var3, :var1), (:var1, :var2, :var3, :var4, :var5, :var6), (:var1, :var2, :var3, :var4, :var5, :var6, :var7, :var8, :var9)],
+            [(:var1,), (:var1, :var2), (:var1, :var3), (:var1, :var4, :var6, :var5), #=(:var1, :var1),=#
+                (:var2, :var7, :var3, :var1), (:var1, :var2, :var3, :var4, :var5, :var6), (:var1, :var2, :var3, :var4, :var5, :var6, :var7, :var0)], #=:var8, :var9,=#
+        ]
 
 
 
@@ -281,17 +276,17 @@ end=#
         for j in 1:length(meteos)
             meteo = meteos[j]
             for k in 1:length(outs_vector)
-                out_tuple  = outs_vector[k]
+                out_tuple = outs_vector[k]
                 #print(i, " ", j, " ", k)
                 meteo_adjusted = PlantSimEngine.adjust_weather_timesteps_to_given_length(
-                    PlantSimEngine.get_status_vector_max_length(modellist.status) , meteo)                             
+                    PlantSimEngine.get_status_vector_max_length(modellist.status), meteo)
                 mtg, mapping, outputs_mapping, nsteps, filtered_outputs_modellist = test_filtered_output_begin(modellist, status_tuple, out_tuple, meteo_adjusted)
                 @test to_initialize(mapping) == Dict()
                 @test test_filtered_output(mtg, mapping, nsteps, outputs_mapping, meteo_adjusted, filtered_outputs_modellist)
             end
         end
     end
-    
+
     #mtg, mapping, outputs_mapping, nsteps, filtered_outputs_modellist = test_filtered_output_begin(modellists[1], status_tuples[1], outs_vectors[1][1], meteos[1])
     #@test test_filtered_output(mtg, mapping, nsteps, outputs_mapping, meteo_day, filtered_outputs_modellist)
 end
@@ -305,7 +300,7 @@ end
 
 function PlantSimEngine.run!(::Reeb, models, status, meteo, constants, extra=nothing)
     status.LAI =
-        status.aPPFD + 0.4*k
+        status.aPPFD + 0.4 * k
 end
 
 function PlantSimEngine.inputs_(::Reeb)
diff --git a/test/test-mapping.jl b/test/test-mapping.jl
index a96c6aa80..96dda7c33 100755
--- a/test/test-mapping.jl
+++ b/test/test-mapping.jl
@@ -72,29 +72,29 @@ end
     mapping_with_vector = Dict(
         "Scale" =>
             (ToyAssimGrowthModel(0.0, 0.0, 0.0),
-            ToyCAllocationModel(),
-            Status( TT_cu=Vector(cumsum(meteo_day.TT))),
+                ToyCAllocationModel(),
+                Status(TT_cu=Vector(cumsum(meteo_day.TT))),
             ),
-        )
-        
-        mtg = import_mtg_example()
-        @test !last(PlantSimEngine.check_statuses_contain_no_remaining_vectors(mapping_with_vector))
-        @test_throws "call the function generate_models_from_status_vectors" PlantSimEngine.GraphSimulation(mtg, mapping_with_vector)
-    
-     mapping_with_empty_status = Dict(
+    )
+
+    mtg = import_mtg_example()
+    @test !last(PlantSimEngine.check_statuses_contain_no_remaining_vectors(mapping_with_vector))
+    @test_throws "call the function generate_models_from_status_vectors" PlantSimEngine.GraphSimulation(mtg, mapping_with_vector)
+
+    mapping_with_empty_status = Dict(
         "Scale" =>
             (ToyAssimGrowthModel(0.0, 0.0, 0.0),
-            ToyCAllocationModel(),
-            Status(),
+                ToyCAllocationModel(),
+                Status(),
             ),
-        )
-    
-     @test last(PlantSimEngine.check_statuses_contain_no_remaining_vectors(mapping_with_empty_status))
+    )
+
+    @test last(PlantSimEngine.check_statuses_contain_no_remaining_vectors(mapping_with_empty_status))
 end
 
 # simple conversion to a mapping, with a manually written model
 function modellist_to_mapping_manual(modellist_original::ModelList, modellist_status, nsteps; check=true, outputs=nothing, TT_cu_vec=Vector{Float64}())
-    
+
     modellist = Base.copy(modellist_original, modellist_original.status)
 
     default_scale = "Default"
@@ -106,17 +106,17 @@ function modellist_to_mapping_manual(modellist_original::ModelList, modellist_st
 
     mapping = Dict(
         default_scale => (
-        models..., 
-        ToyTestDegreeDaysCumulModel(TT_cu_vec=TT_cu_vec),
-        PlantSimEngine.HelperNextTimestepModel(),
-        MultiScaleModel(
-        model=PlantSimEngine.HelperCurrentTimestepModel(),
-        mapped_variables=[PreviousTimeStep(:next_timestep),],
-        ),
-        Status(current_timestep=1,next_timestep=1)
+            models...,
+            ToyTestDegreeDaysCumulModel(TT_cu_vec=TT_cu_vec),
+            PlantSimEngine.HelperNextTimestepModel(),
+            MultiScaleModel(
+                model=PlantSimEngine.HelperCurrentTimestepModel(),
+                mapped_variables=[PreviousTimeStep(:next_timestep),],
+            ),
+            Status(current_timestep=1, next_timestep=1)
         ),
     )
-    
+
     if isnothing(outputs)
         f = []
         for i in 1:length(modellist.models)
@@ -131,7 +131,7 @@ function modellist_to_mapping_manual(modellist_original::ModelList, modellist_st
         f = unique!(f)
         all_vars = (f...,)
         #all_vars = merge((keys(init_variables(object.models[i])) for i in 1:length(object.models))...)
-    else 
+    else
         all_vars = outputs
         # TODO sanity check
     end
@@ -150,7 +150,7 @@ end
 PlantSimEngine.inputs_(::ToyTestDegreeDaysCumulModel) = (current_timestep=1,)
 PlantSimEngine.outputs_(::ToyTestDegreeDaysCumulModel) = (TT_cu=0.0,)
 
-ToyTestDegreeDaysCumulModel(; TT_cu_vec = Vector{Float64}()) = ToyTestDegreeDaysCumulModel(TT_cu_vec)
+ToyTestDegreeDaysCumulModel(; TT_cu_vec=Vector{Float64}()) = ToyTestDegreeDaysCumulModel(TT_cu_vec)
 
 
 function PlantSimEngine.run!(m::ToyTestDegreeDaysCumulModel, models, status, meteo, constants=nothing, extra=nothing)
@@ -164,7 +164,7 @@ PlantSimEngine.ObjectDependencyTrait(::Type{<:ToyTestDegreeDaysCumulModel}) = Pl
     meteo_day = CSV.read(joinpath(pkgdir(PlantSimEngine), "examples/meteo_day.csv"), DataFrame, header=18)
 
     st = (TT_cu=cumsum(meteo_day.TT),)
-    
+
     TT_cu_vec = Vector(cumsum(meteo_day.TT))
 
     rue = 0.3
@@ -197,9 +197,9 @@ PlantSimEngine.ObjectDependencyTrait(::Type{<:ToyTestDegreeDaysCumulModel}) = Pl
 
     # fully automated model generation
     st2 = (TT_cu=Vector(cumsum(meteo_day.TT)),)
-   
+
     mtg, mapping, outputs_mapping = PlantSimEngine.modellist_to_mapping(models, st2; nsteps=nsteps, outputs=nothing)
- 
+
     @test to_initialize(mapping) == Dict()
 
     graphsim2 = PlantSimEngine.GraphSimulation(mtg, mapping, nsteps=nsteps, check=true, outputs=outputs_mapping)
@@ -225,52 +225,50 @@ end
     TT_cu_vec = Vector(cumsum(meteo_day.TT))
     nsteps = length(meteo_day.TT)
 
-    mapping_with_vector = Dict(
-    
-    "Plant" => (
-        MultiScaleModel(
-            model=ToyCAllocationModel(),
-            mapped_variables=[
-                # inputs
-                :carbon_assimilation => ["Leaf"],
-                :carbon_demand => ["Leaf", "Internode"],
-                # outputs
-                :carbon_allocation => ["Leaf", "Internode"]
-            ],
+    mapping_with_vector = Dict("Plant" => (
+            MultiScaleModel(
+                model=ToyCAllocationModel(),
+                mapped_variables=[
+                    # inputs
+                    :carbon_assimilation => ["Leaf"],
+                    :carbon_demand => ["Leaf", "Internode"],
+                    # outputs
+                    :carbon_allocation => ["Leaf", "Internode"]
+                ],
+            ),
+            MultiScaleModel(
+                model=ToyPlantRmModel(),
+                mapped_variables=[:Rm_organs => ["Leaf" => :Rm, "Internode" => :Rm],],
+            ),
         ),
-        MultiScaleModel(
-            model=ToyPlantRmModel(),
-            mapped_variables=[:Rm_organs => ["Leaf" => :Rm, "Internode" => :Rm],],
+        "Internode" => (
+            ToyCDemandModel(optimal_biomass=10.0, development_duration=200.0),
+            ToyMaintenanceRespirationModel(1.5, 0.06, 25.0, 0.6, 0.004),
+            Status(TT=TT_v, carbon_biomass=1.0)
         ),
-    ),
-    "Internode" => (
-        ToyCDemandModel(optimal_biomass=10.0, development_duration=200.0),
-        ToyMaintenanceRespirationModel(1.5, 0.06, 25.0, 0.6, 0.004),
-        Status(TT=TT_v, carbon_biomass=1.0)
-    ),
-    "Leaf" => (
-        MultiScaleModel(
-            model=ToyAssimModel(),
-            mapped_variables=[:soil_water_content => "Soil",],
-            # Notice we provide "Soil", not ["Soil"], so a single value is expected here
+        "Leaf" => (
+            MultiScaleModel(
+                model=ToyAssimModel(),
+                mapped_variables=[:soil_water_content => "Soil",],
+                # Notice we provide "Soil", not ["Soil"], so a single value is expected here
+            ),
+            ToyCDemandModel(optimal_biomass=10.0, development_duration=200.0),
+            ToyMaintenanceRespirationModel(2.1, 0.06, 25.0, 1.0, 0.025),
+            Status(aPPFD=1300.0, carbon_biomass=2.0, TT=10.0), # TODO try calling the generated TT output through a variable mapping
         ),
-        ToyCDemandModel(optimal_biomass=10.0, development_duration=200.0),
-        ToyMaintenanceRespirationModel(2.1, 0.06, 25.0, 1.0, 0.025),
-        Status(aPPFD=1300.0, carbon_biomass=2.0, TT=10.0), # TODO try calling the generated TT output through a variable mapping
-    ),
-    "Soil" => (
-        ToySoilWaterModel(),
-    ),
-)
+        "Soil" => (
+            ToySoilWaterModel(),
+        ),
+    )
 
-out_multiscale = Dict("Plant" => (:Rm_organs,),)
-mtg = import_mtg_example();
+    out_multiscale = Dict("Plant" => (:Rm_organs,),)
+    mtg = import_mtg_example()
 
-mapping_without_vectors = PlantSimEngine.replace_mapping_status_vectors_with_generated_models(mapping_with_vector, "Soil", nsteps)
+    mapping_without_vectors = PlantSimEngine.replace_mapping_status_vectors_with_generated_models(mapping_with_vector, "Soil", nsteps)
 
-@test to_initialize(mapping_without_vectors) == Dict()
+    @test to_initialize(mapping_without_vectors) == Dict()
 
- graph_sim_multiscale = @test_nowarn PlantSimEngine.GraphSimulation(mtg, mapping_without_vectors, nsteps=nsteps, check=true, outputs=out_multiscale)
+    graph_sim_multiscale = @test_nowarn PlantSimEngine.GraphSimulation(mtg, mapping_without_vectors, nsteps=nsteps, check=true, outputs=out_multiscale)
 
     sim_multiscale = run!(graph_sim_multiscale,
         meteo_day,