Repository: SDXorg/pysd (language: python)

Issue:
get_dependencies doc does not match its behaviour
In the description of the Model.get_dependencies documentation it says that the method will return a set, but actually there is no return, only dependencies printed out.
https://github.com/SDXorg/pysd/blob/bcfa21141c013c11028c611be1e627411c6b85d3/pysd/py_backend/model.py#L1376 
What would be the expected behaviour for that method ?

Only printing or actually returning the dependencies ?

I would suggest to simply return dependencies as a dependencies object and have the __str__ method being what is currently printed and supporting iteration and accessing similarly to a python set.
Implementing externals serialization for faster model initialization
## Description

This PR adds the possibility to serialize external data in the form of a netCDF file, and then use the file to load external data from it.

Also adds importlib_metadata in requirements.txt because xarray imports entry_points() from importlib_metadata instead of importlib.metadata for Python < 3.8.

## Related issues

Closes #355 

## Type of change

- New feature (non-breaking change which adds functionality)

## PR verification (to be filled by reviewers)

- [ ] The code follows the [PEP 8 style](https://peps.python.org/pep-0008/)
- [ ] The new code has been tested properly for all the possible cases
- [ ] The overall coverage has not dropped and other features have not been broken.
- [ ] If new features have been added, they have been properly documented
- [ ] *docs/whats_new.rst* has been updated
- [ ] PySD version has been updated

Add Dependabot configuration for GitHub Actions updates
Add a Dependabot configuration that checks once a week if the GitHub Actions are still using the latest version. If not, it opens a PR to update them.

It will actually open very few PRs, since we only have major versions specified (like v3), so only on a major v4 release it will update and open a PR.

See [Keeping your actions up to date with Dependabot](https://docs.github.com/en/code-security/dependabot/working-with-dependabot/keeping-your-actions-up-to-date-with-dependabot).
Repeated x values in HardcodedLookups produces "InvalidIndexError: Reindexing only valid with uniquely valued Index objects after many loops"
@ali-akhavan89 shared the following error and we think this maybe because of the structure below (when `inventory coverage` is removed, Ali said the error disappered), but am not sure.

<img width="712" alt="image" src="https://user-images.githubusercontent.com/30194633/194606692-2ec19aa5-13fe-450a-b210-503abb56a075.png">


I think For the attached `.mdl` file, 
> model = pysd.read_vensim('test.mdl')
> values = model.run()

text for `test.mdl`
```

Desired Shipment=
	100
	~	Units/Year
	~		|

Inventory= INTEG (
	Production Rate-Shipment Rate,
		20)
	~	Units
	~		|

Inventory Coverage=
	Inventory/Shipment Rate
	~	Year
	~		|

Maximum Shipment Rate=
	Inventory/Minimum Order Processing Time
	~	Units/Years
	~		|

Minimum Order Processing Time=
	0.1
	~	Years
	~		|

Order Fulfillment Ratio=
	Table for Order Fulfillment(Maximum Shipment Rate/Desired Shipment)
	~	Dmnl
	~		|

Production Rate=
	100
	~	Units/Year
	~		|

Shipment Rate=
	Desired Shipment*Order Fulfillment Ratio
	~	Units/Year
	~		|

Table for Order Fulfillment(
	[(0,0)-(4,1)],(0,0),(0.2,0.2),(0.4,0.4),(0.6,0.58),(0.8,0.73),(1,0.85),(1.2,0.93),(1.4\
		,0.97),(1.6,0.99),(1.8,1),(2,1),(2,1),(4,1))
	~	Dmnl
	~		|

********************************************************
	.Control
********************************************************~
		Simulation Control Parameters
	|

FINAL TIME  = 160
	~	Year
	~	The final time for the simulation.
	|

INITIAL TIME  = 0
	~	Year
	~	The initial time for the simulation.
	|

SAVEPER  = 
        TIME STEP
	~	Year [0,?]
	~	The frequency with which output is stored.
	|

TIME STEP  = 0.03125
	~	Year [0,?]
	~	The time step for the simulation.
	|

```
Error log:
```
<locals>.cached_func(*args)
24 try:
---> 25 return self.data[func.__name__]
26 except KeyError:
KeyError: 'shipment_rate'
During handling of the above exception, another exception
occurred:
KeyError Traceback (most
recent call last)

File ~/opt/anaconda3/lib/python3.9/site-
packages/pysd/py_backend/cache.py:25, in Cache.__call__.

<locals>.cached_func(*args)
24 try:
---> 25 return self.data[func.__name__]
26 except KeyError:
KeyError: 'order_fulfillment_ratio'
During handling of the above exception, another exception
occurred:
InvalidIndexError Traceback (most
recent call last)
Input In [1], in <cell line: 20>()
16 import xarray
19 model =
pysd.read_vensim('/Users/Ali/PySD/Commodity/test10.mdl')
---> 20 values = model.run()
21 values.head(10)

File ~/opt/anaconda3/lib/python3.9/site-
packages/pysd/py_backend/model.py:1214, in Model.run(self,

params, return_columns, return_timestamps,
initial_condition, final_time, time_step, saveper, reload,
progress, flatten_output, cache_output, output_file)
1211 # instantiating output object

1212 output = ModelOutput(self, capture_elements['step'],
output_file)
-> 1214 self._integrate(output)
1216 del self._dependencies["OUTPUTS"]
1218 output.add_run_elements(self,
capture_elements['run'])

File ~/opt/anaconda3/lib/python3.9/site-
packages/pysd/py_backend/model.py:1719, in

Model._integrate(self, out_obj)
1717 while self.time.in_bounds():
1718 if self.time.in_return():
-> 1719 out_obj.update(self)
1721 self._euler_step(self.time.time_step())
1722
self.time.update(self.time()+self.time.time_step())

File ~/opt/anaconda3/lib/python3.9/site-
packages/pysd/py_backend/output.py:59, in

ModelOutput.update(self, model)
56 def update(self, model):
57 """ Delegating the update of the results object
and its elements to the
58 appropriate handler."""
---> 59 self.handler.update(model,
self.capture_elements)

File ~/opt/anaconda3/lib/python3.9/site-
packages/pysd/py_backend/output.py:442, in

DataFrameHandler.update(self, model, capture_elements)
425 def update(self, model, capture_elements):
426 """
427 Add a row to the results pandas DataFrame with the
values of the
428 variables listed in capture_elements.
(...)
440
441 """
--> 442 self.ds.loc[model.time.round()] = [
443 getattr(model.components, key)()
444 for key in capture_elements]

File ~/opt/anaconda3/lib/python3.9/site-
packages/pysd/py_backend/output.py:443, in <listcomp>(.0)

425 def update(self, model, capture_elements):

426 """
427 Add a row to the results pandas DataFrame with the
values of the
428 variables listed in capture_elements.
(...)
440
441 """
442 self.ds.loc[model.time.round()] = [
--> 443 getattr(model.components, key)()
444 for key in capture_elements]
File ~/PySD/Commodity/test10.py:141, in
inventory_coverage()
133 @component.add(
134 name="Inventory Coverage",
135 units="Year",
(...)
139 )
140 def inventory_coverage():
--> 141 return inventory() / shipment_rate()

File ~/opt/anaconda3/lib/python3.9/site-
packages/pysd/py_backend/cache.py:27, in Cache.__call__.

<locals>.cached_func(*args)
25 return self.data[func.__name__]
26 except KeyError:
---> 27 value = func(*args)
28 self.data[func.__name__] = value
29 return value
File ~/PySD/Commodity/test10.py:198, in shipment_rate()
190 @component.add(
191 name="Shipment Rate",
192 units="Units/Year",
(...)
196 )
197 def shipment_rate():
--> 198 return desired_shipment() *
order_fulfillment_ratio()

File ~/opt/anaconda3/lib/python3.9/site-
packages/pysd/py_backend/cache.py:27, in Cache.__call__.

<locals>.cached_func(*args)
25 return self.data[func.__name__]
26 except KeyError:

---> 27 value = func(*args)
28 self.data[func.__name__] = value
29 return value
File ~/PySD/Commodity/test10.py:177, in
order_fulfillment_ratio()
165 @component.add(
166 name="Order Fulfillment Ratio",
167 units="Dmnl",
(...)
175 )
176 def order_fulfillment_ratio():
--> 177 return
table_for_order_fulfillment(maximum_shipment_rate() /
desired_shipment())
File ~/PySD/Commodity/test10.py:209, in
table_for_order_fulfillment(x, final_subs)
201 @component.add(
202 name="Table for Order Fulfillment",
203 units="Dmnl",
(...)
207 )
208 def table_for_order_fulfillment(x, final_subs=None):
--> 209 return
_hardcodedlookup_table_for_order_fulfillment(x, final_subs)

File ~/opt/anaconda3/lib/python3.9/site-
packages/pysd/py_backend/lookups.py:39, in

Lookups.__call__(self, x, final_subs)
37 def __call__(self, x, final_subs=None):
38 try:
---> 39 return self._call(self.data, x, final_subs)
40 except (TypeError, KeyError):
41 # this except catch the errors when a lookups
has been
42 # changed to a constant value by the user
43 if final_subs and isinstance(self.data,
xr.DataArray):
44 # self.data is an array, reshape it

File ~/opt/anaconda3/lib/python3.9/site-
packages/pysd/py_backend/lookups.py:87, in

Lookups._call(self, data, x, final_subs)
85 else:

86 if x in data['lookup_dim'].values:
---> 87 outdata = data.sel(lookup_dim=x)
88 elif x > data['lookup_dim'].values[-1]:
89 if self.interp == "extrapolate":
90 # extrapolate method for xmile models

File ~/opt/anaconda3/lib/python3.9/site-
packages/xarray/core/dataarray.py:1332, in

DataArray.sel(self, indexers, method, tolerance, drop,
**indexers_kwargs)
1223 def sel(
1224 self,
1225 indexers: Mapping[Any, Any] = None,
(...)
1229 **indexers_kwargs: Any,
1230 ) -> "DataArray":
1231 """Return a new DataArray whose data is given by
selecting index
1232 labels along the specified dimension(s).
1233
(...)
1330 Dimensions without coordinates: points
1331 """
-> 1332 ds = self._to_temp_dataset().sel(
1333 indexers=indexers,
1334 drop=drop,
1335 method=method,
1336 tolerance=tolerance,
1337 **indexers_kwargs,
1338 )
1339 return self._from_temp_dataset(ds)

File ~/opt/anaconda3/lib/python3.9/site-
packages/xarray/core/dataset.py:2504, in Dataset.sel(self,

indexers, method, tolerance, drop, **indexers_kwargs)
2443 """Returns a new dataset with each array indexed by
tick labels
2444 along the specified dimension(s).
2445
(...)
2501 DataArray.sel
2502 """
2503 indexers = either_dict_or_kwargs(indexers,
indexers_kwargs, "sel")
-> 2504 pos_indexers, new_indexes = remap_label_indexers(

2505 self, indexers=indexers, method=method,
tolerance=tolerance
2506 )
2507 # TODO: benbovy - flexible indexes: also use
variables returned by Index.query
2508 # (temporary dirty fix).
2509 new_indexes = {k: v[0] for k, v in new_indexes.items()}

File ~/opt/anaconda3/lib/python3.9/site-
packages/xarray/core/coordinates.py:421, in

remap_label_indexers(obj, indexers, method, tolerance,
**indexers_kwargs)
414 indexers = either_dict_or_kwargs(indexers,
indexers_kwargs, "remap_label_indexers")
416 v_indexers = {
417 k: v.variable.data if isinstance(v, DataArray) else
v
418 for k, v in indexers.items()
419 }
--> 421 pos_indexers, new_indexes =
indexing.remap_label_indexers(
422 obj, v_indexers, method=method,
tolerance=tolerance
423 )
424 # attach indexer's coordinate to pos_indexers
425 for k, v in indexers.items():

File ~/opt/anaconda3/lib/python3.9/site-
packages/xarray/core/indexing.py:120, in

remap_label_indexers(data_obj, indexers, method, tolerance)
118 for dim, index in indexes.items():
119 labels = grouped_indexers[dim]
--> 120 idxr, new_idx = index.query(labels,
method=method, tolerance=tolerance)
121 pos_indexers[dim] = idxr
122 if new_idx is not None:

File ~/opt/anaconda3/lib/python3.9/site-
packages/xarray/core/indexes.py:240, in

PandasIndex.query(self, labels, method, tolerance)
238 indexer = label
239 else:
--> 240 indexer = get_indexer_nd(self.index, label,
method, tolerance)
241 if np.any(indexer < 0):

242 raise KeyError(f"not all values found in index
{coord_name!r}")

File ~/opt/anaconda3/lib/python3.9/site-
packages/xarray/core/indexes.py:142, in

get_indexer_nd(index, labels, method, tolerance)
138 """Wrapper around :meth:`pandas.Index.get_indexer`
supporting n-dimensional
139 labels
140 """
141 flat_labels = np.ravel(labels)
--> 142 flat_indexer = index.get_indexer(flat_labels,
method=method, tolerance=tolerance)
143 indexer = flat_indexer.reshape(labels.shape)
144 return indexer

File ~/opt/anaconda3/lib/python3.9/site-
packages/pandas/core/indexes/base.py:3721, in

Index.get_indexer(self, target, method, limit, tolerance)
3718 self._check_indexing_method(method, limit, tolerance)
3720 if not self._index_as_unique:
-> 3721 raise
InvalidIndexError(self._requires_unique_msg)
3723 if len(target) == 0:
3724 return np.array([], dtype=np.intp)
InvalidIndexError: Reindexing only valid with uniquely
valued Index objects
```
get_dependencies doc does not match its behaviour
In the description of the Model.get_dependencies documentation it says that the method will return a set, but actually there is no return, only dependencies printed out.
https://github.com/SDXorg/pysd/blob/bcfa21141c013c11028c611be1e627411c6b85d3/pysd/py_backend/model.py#L1376 
What would be the expected behaviour for that method ?

Only printing or actually returning the dependencies ?

I would suggest to simply return dependencies as a dependencies object and have the __str__ method being what is currently printed and supporting iteration and accessing similarly to a python set.
Allow exporting and reading netCDF data
Some models may require to read large quantities of data from spreadsheets. This may  provoke a slow model initialization. When this data is not an user expected input but some constants that we are not expecting to change from one simulation to another, it would be feasible to move it to a binary file. As we are supporting multi dimensional arrays with xarray, I thik that the most appropriate file format to use would be netCDF. xarray allows easily reading and exporting to netCDF, moreover there are some other tools to read this data. We should:

- Provide methods to export all or some External objects data to netCDF.
- Provide methods to read External objects data  from netCDF, taking into account that maybe some External objects should still be read from spreadsheets.
- Think if we want to have the expression for reading the data from a excel file or from a netCDF file in the final model. The second one will produce mutch simpler code, but require to have always the netCDF files available.


Discussion / hints:

Fix patch:
diff --git a/.github/dependabot.yml b/.github/dependabot.yml
new file mode 100644
index 00000000..ac27a848
--- /dev/null
+++ b/.github/dependabot.yml
@@ -0,0 +1,8 @@
+version: 2
+updates:
+  # Maintain dependencies for GitHub Actions
+  - package-ecosystem: "github-actions"
+    directory: "/"
+    schedule:
+      # Check for updates to GitHub Actions every week
+      interval: "weekly"
diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 6d51229e..f9981634 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -25,7 +25,7 @@ jobs:
           pip install -r tests/requirements.txt
           pip install -e .
       - name: Test and coverage
-        run: pytest tests/ --cov=pysd
+        run: pytest tests/ --cov=pysd -n 2
       - name: Coveralls
         if: ${{ matrix.python-version == 3.7 && matrix.os == 'ubuntu-latest' }}
         env:
diff --git a/docs/advanced_usage.rst b/docs/advanced_usage.rst
index 9232c35c..debe650d 100644
--- a/docs/advanced_usage.rst
+++ b/docs/advanced_usage.rst
@@ -125,3 +125,63 @@ In order to preview the needed exogenous variables, the :py:meth:`.get_dependenc
 
 .. automethod:: pysd.py_backend.model.Model.get_dependencies
    :noindex:
+
+
+Initializing external data from netCDF file
+-------------------------------------------
+
+IO operations are expensive, especially when reading non-binary files. This makes the model initialization process slow when lots of datasets need to be read from spreadhseet files.
+
+From PySD 3.8, users can export a subset or all model external data to a netCDF file, and use this file for subsequent model initializations.
+
+Supose we have a model (*my_model.mdl*) that loads *param_1* from *parameters_1.xls*, *param_2* from *parameters_2.xls*, and *policy_1* and *policy_2* from *scenario.xls*. Imagine we want to test different policy configurations, by changing the values of *policy_1* and *policy_2*, while keeping all other parameters unchanged.
+In this case, we might want to export the external objects that we do not intend to modify (*param_1* and *param_2*) to a netCDF file, so that they are initialized instantaneously:
+
+The code below shows how this can be achieved::
+
+    import pysd
+    model = pysd.read_vensim("my_model.mdl", initialize=False)
+    model.serialize_externals(export_path="parameters.nc",
+                              include_externals="all",
+                              exclude_externals=["scenario.xls"])
+
+
+This will export only the external data defined in *parameters_1.xls* and *parameters_2.xls* to the *parameters.xls* file. A good thing of this is that the resulting netCDF file will include all the metadata defined in the model for such variable (description, units, etc.)
+
+Note that the exact same result could have been achieved with::
+
+    import pysd
+
+    model = pysd.read_vensim("my_model.mdl", initialize=False)
+    model.serialize_externals(export_path="parameters.nc",
+                              include_externals="all",
+                              exclude_externals=["policy_1", "policy_2"])
+
+Or even::
+
+    import pysd
+
+    model = pysd.read_vensim("my_model.mdl", initialize=False)
+    model.serialize_externals(export_path="parameters.nc",
+                              include_externals=["parameters_1.xls",
+                                                 "parameters_2.xls"],
+                              exclude_externals=None)
+
+Or we could have also combined variable names and spreadhseet files in the **include_externals** arugment, the **exclude_externals** argument or both::
+
+    import pysd
+
+    model = pysd.read_vensim("my_model.mdl", initialize=False)
+    model.serialize_externals(export_path="parameters.nc",
+                              include_externals=["param_1", "parameters_2.xls"],
+                              exclude_externals=None)
+
+
+Then, to run a simulation loading the external objects stored in *parameters.nc*, we write::
+
+    import pysd
+
+    model = pysd.read_vensim("my_model.mdl", initialize=False)
+    model.initialize_external_data(externals="parameters.nc")
+
+And from here we could run the model normally with :py:meth:`pysd.py_backend.Model.run`.
\ No newline at end of file
diff --git a/docs/getting_started.rst b/docs/getting_started.rst
index bfb739bc..e4257866 100644
--- a/docs/getting_started.rst
+++ b/docs/getting_started.rst
@@ -175,8 +175,11 @@ If the `output_file` is not set, the :py:meth:`.run` method will return a :py:cl
 
 For most cases, the *.tab* file format is the safest choice. It is preferable over the *.csv* format when the model includes subscripted variables. The *.nc* format is recommended for large models, and when the user wants to keep metadata such as variable units and description.
 
+.. note::
+   PySD includes helper functions to export *.nc* file contents to *.csv* or *.tab* files. See :doc:`Exporting netCDF data_vars to csv or tab <./tools>` for further details.
+
 .. warning::
-   *.nc* files require :py:mod:`netcdf4` library which is an optional requirement and thus not installed automatically with the package. We recommend using :py:mod:`netcdf4` 1.6.0 or above, however, it will also work with :py:mod:`netcdf4` 1.5.0 or above.
+   *.nc* files require :py:mod:`netCDF4` library which is an optional requirement and thus not installed automatically with the package. We recommend using :py:mod:`netCDF4` 1.6.0 or above, however, it will also work with :py:mod:`netCDF4` 1.5.0 or above.
 
 
 Setting parameter values
diff --git a/docs/installation.rst b/docs/installation.rst
index 6f5cf8ed..6ba9a680 100644
--- a/docs/installation.rst
+++ b/docs/installation.rst
@@ -66,9 +66,16 @@ In order to plot model outputs as shown in :doc:`Getting started <../getting_sta
 
 * Matplotlib
 
-In order to be able to export data to netCDF (*.nc*) files:
+To export data to netCDF (*.nc*) files:
 
 * netCDF4
+* importlib-metadata==2.0 (only for Python 3.7)
+
+To export netCDF data to comma or tab separated files with parallel processing:
+
+* dask[array]
+* dask[diagnostics]
+* dask[distributed]
 
 
 These Python libraries bring additional data analytics capabilities to the analysis of SD models:
diff --git a/docs/tools.rst b/docs/tools.rst
index 99aa61f5..d978c191 100644
--- a/docs/tools.rst
+++ b/docs/tools.rst
@@ -7,3 +7,20 @@ Benchmarking
 ------------
 .. automodule:: pysd.tools.benchmarking
    :members:
+
+Exporting netCDF data_vars to csv or tab
+----------------------------------------
+
+Simulation results can be stored as netCDF (*.nc*) files (see :doc:`Storing simulation results on a file <./getting_started>`).
+
+The :py:class:`pysd.tools.ncfiles.NCFile` allows loading netCDF files generated with PySD as an :py:class:`xarray.Dataset`. When passing the argument `parallel=True` to the constructor, :py:class:`xarray.DataArray` inside the Dataset will be loded as `dask arrays <https://docs.dask.org/en/stable/array.html>`_, with `chunks=-1 <https://docs.dask.org/en/stable/array-chunks.html>`_.
+
+Once the Dataset is loaded, a subset (or all) of the data_vars can be exported into:
+
+* A :py:class:`pandas.DataFrame`, using the :py:meth:`pysd.tools.ncfiles.NCFile.to_df` method
+* A `*.csv` or `*.tab` files, using the the :py:meth:`pysd.tools.ncfiles.NCFile.to_text_file` method
+
+Alternatively, to get further control of the chunking, users can load the :py:class:`xarray.Dataset` using :py:meth:`xarray.open_dataset` and then use the :py:meth:`pysd.tools.ncfiles.NCFile.ds_to_df` or :py:meth:`pysd.tools.ncfiles.NCFile.df_to_text_file` static methods.
+
+.. automodule:: pysd.tools.ncfiles
+   :members:
\ No newline at end of file
diff --git a/docs/whats_new.rst b/docs/whats_new.rst
index 035f498d..5e4719a6 100644
--- a/docs/whats_new.rst
+++ b/docs/whats_new.rst
@@ -1,5 +1,44 @@
 What's New
 ==========
+v3.8.0 (2022/11/03)
+-------------------
+
+New Features
+~~~~~~~~~~~~
+- Adds ncfile.py module with helper functions to export a subset or all of the data_vars in netCDF files generated with PySD to :py:class:`pandas.DataFrame`, csv or tab files. (`@rogersamso <https://github.com/rogersamso>`_)
+- Adds possibility to initialize and export a subset or all external objects to netCDF, and then initialize the external objects from the file. (`@rogersamso <https://github.com/rogersamso>`_)
+
+Breaking changes
+~~~~~~~~~~~~~~~~
+
+Deprecations
+~~~~~~~~~~~~
+- Deprecate :py:meth:`pysd.py_backend.model.Model._get_dependencies` replacing it with :py:meth:`pysd.py_backend.model.Model.get_dependencies`. (`@enekomartinmartinez <https://github.com/enekomartinmartinez>`_)
+
+Bug fixes
+~~~~~~~~~
+- Include new class :py:class:`pysd.py_backend.utils.Dependencies` to return by :py:meth:`pysd.py_backend.model.Model.get_dependencies` (:issue:`379`). (`@lionel42 <https://github.com/lionel42>`_)
+
+Documentation
+~~~~~~~~~~~~~
+- Updates the :doc:`getting_started` page with instructions on how to use the new helper functions for netCDF files. (`@rogersamso <https://github.com/rogersamso>`_)
+- Updates the :doc:`advanced_usage` page with instructions on how to export externals to netCDF and initialize a model from it. (`@rogersamso <https://github.com/rogersamso>`_)
+- Update citation information to include the new paper published in JOSS. (`@enekomartinmartinez <https://github.com/enekomartinmartinez>`_)
+
+Performance
+~~~~~~~~~~~
+- Initializing external objects from netCDF file is much faster than reading from spreadsheet files.(`@rogersamso <https://github.com/rogersamso>`_)
+
+Internal Changes
+~~~~~~~~~~~~~~~~
+- Adds the :py:meth:`pysd.py_backend.model.Macro.serialize_externals` and :py:meth:`pysd.py_backend.model.Macro.initialize_external_data` methods, and a few other private methods.(`@rogersamso <https://github.com/rogersamso>`_)
+- Adds the :py:class:`pysd.py_backend.utils.UniqueDims` class for renaming model dimensions with unique names.(`@rogersamso <https://github.com/rogersamso>`_)
+- Force :py:class:`pysd.py_backend.external.External` objects to always have the full element dimensions, missing dimensions are filled with `numpy.nan`. (`@enekomartinmartinez <https://github.com/enekomartinmartinez>`_)
+- Add `dependabot <https://docs.github.com/en/code-security/dependabot/working-with-dependabot/keeping-your-actions-up-to-date-with-dependabot>`_ configuration for GitHub Actions updates. (`@EwoutH <https://github.com/EwoutH>`_)
+- Include new error messages for initialization of :py:class:`pysd.py_backend.lookups.HardcodedLookups` (:issue:`376`). (`@enekomartinmartinez <https://github.com/enekomartinmartinez>`_)
+- Include new warning message when a translated variable has several types or subtypes. (`@enekomartinmartinez <https://github.com/enekomartinmartinez>`_)
+- Set CI test to run in parallel in 2 cores. (`@enekomartinmartinez <https://github.com/enekomartinmartinez>`_)
+
 v3.7.1 (2022/09/19)
 -------------------
 
diff --git a/pysd/_version.py b/pysd/_version.py
index 975f6914..32a78190 100644
--- a/pysd/_version.py
+++ b/pysd/_version.py
@@ -1,1 +1,1 @@
-__version__ = "3.7.1"
+__version__ = "3.8.0"
diff --git a/pysd/builders/python/python_expressions_builder.py b/pysd/builders/python/python_expressions_builder.py
index e4df5a51..753e2388 100644
--- a/pysd/builders/python/python_expressions_builder.py
+++ b/pysd/builders/python/python_expressions_builder.py
@@ -279,7 +279,8 @@ def update_object_subscripts(self, name: str,
                                  component_final_subs: dict) -> None:
         """
         Update the object subscripts. Needed for those objects that
-        use 'add' method to load several components at once.
+        use 'add' method to load several components at once and mixed
+        definitions are used.
 
         Parameters
         ----------
@@ -295,39 +296,8 @@ def update_object_subscripts(self, name: str,
         """
         # Get the component used to define the object first time
         origin_comp = self.element.objects[name]["component"]
-        if isinstance(origin_comp.subscripts_dict, dict):
-            # The original component subscript dictionary is a dict
-            if len(list(origin_comp.subscripts_dict)) == 1:
-                # If the subscript dict has only one dimension
-                # all the components can be loaded in 1D array directly
-                # with the same length as the given by the sum of the
-                # components
-                key = list(origin_comp.subscripts_dict.keys())[0]
-                value = list(component_final_subs.values())[0]
-                origin_comp.subscripts_dict[key] += value
-                # force the same order as the final element
-                origin_comp.subscripts_dict[key] = [
-                    sub for sub in self.element.subs_dict[key]
-                    if sub in origin_comp.subscripts_dict[key]
-                ]
-                self.element.objects[name]["final_subs"] =\
-                    origin_comp.subscripts_dict
-            else:
-                # If the subscripts dict has several dimensions, then
-                # a multi-dimensional array needs to be computed,
-                # in some cases, when mixed definitions are used in an
-                # element (e.g. GET DIRECT CONSTANTS and regular constants),
-                # this array can have some empty subarrays, therefore a
-                # list should be created and manage the loaded data
-                # with manage_multi_def in the element building
-                origin_comp.subscripts_dict = [origin_comp.subscripts_dict]
-                self.element.objects[name]["final_subs"] =\
-                    self.element.subs_dict
-        if isinstance(origin_comp.subscripts_dict, list):
-            # The original component subscript dictionary is a list
-            # (this happens when other components have already been
-            # added with 'add' method)
-            origin_comp.subscripts_dict.append(component_final_subs)
+        # The original component subscript dictionary is a list
+        origin_comp.subscripts_dict.append(component_final_subs)
 
 
 class OperationBuilder(StructureBuilder):
@@ -939,8 +909,10 @@ def build(self, arguments: dict) -> Union[BuildAST, None]:
 
             arguments["name"] = self.section.namespace.make_python_identifier(
                 self.element.identifier, prefix="_ext_lookup")
-            arguments["final_subs"] = "%(final_subs)s"
-            self.component.subscripts_dict = final_subs
+            arguments["final_subs"] =\
+                self.section.subscripts.simplify_subscript_input(
+                    self.element.subs_dict)[1]
+            self.component.subscripts_dict = [final_subs]
 
             self.element.objects["ext_lookups"] = {
                 "name": arguments["name"],
@@ -1013,8 +985,10 @@ def build(self, arguments: dict) -> Union[BuildAST, None]:
 
             arguments["name"] = self.section.namespace.make_python_identifier(
                 self.element.identifier, prefix="_ext_data")
-            arguments["final_subs"] = "%(final_subs)s"
-            self.component.subscripts_dict = final_subs
+            arguments["final_subs"] =\
+                self.section.subscripts.simplify_subscript_input(
+                    self.element.subs_dict)[1]
+            self.component.subscripts_dict = [final_subs]
 
             self.element.objects["ext_data"] = {
                 "name": arguments["name"],
@@ -1085,8 +1059,10 @@ def build(self, arguments: dict) -> Union[BuildAST, None]:
 
             arguments["name"] = self.section.namespace.make_python_identifier(
                 self.element.identifier, prefix="_ext_constant")
-            arguments["final_subs"] = "%(final_subs)s"
-            self.component.subscripts_dict = final_subs
+            arguments["final_subs"] =\
+                self.section.subscripts.simplify_subscript_input(
+                    self.element.subs_dict)[1]
+            self.component.subscripts_dict = [final_subs]
 
             self.element.objects["constants"] = {
                 "name": arguments["name"],
@@ -1712,7 +1688,9 @@ def build(self, arguments: dict) -> Union[BuildAST, None]:
 
             arguments["name"] = self.section.namespace.make_python_identifier(
                 self.element.identifier, prefix="_hardcodedlookup")
-            arguments["final_subs"] = "%(final_subs)s"
+            arguments["final_subs"] =\
+                self.section.subscripts.simplify_subscript_input(
+                    self.element.subs_dict)[1]
 
             self.element.objects["hardcoded_lookups"] = {
                 "name": arguments["name"],
diff --git a/pysd/builders/python/python_model_builder.py b/pysd/builders/python/python_model_builder.py
index 637571f6..ce432835 100644
--- a/pysd/builders/python/python_model_builder.py
+++ b/pysd/builders/python/python_model_builder.py
@@ -638,12 +638,7 @@ def build_element(self) -> None:
             self.expression = "value"
         else:
             self.pre_expression = ""
-            # NUMPY: reshape to the final shape if needed
-            # include full subscript range for objects defined with a
-            # partial range (issue #363)
-            for value in self.objects.values():
-                if value["expression"] is not None and "final_subs" in value:
-                    value["final_subs"] = self.subs_dict
+            # NUMPY: reshape to the final shape if needed)
             # expressions[0]["expr"].reshape(self.section.subscripts, {})
             if not expressions[0]["expr"].subscripts and self.subscripts:
                 # Updimension the return value to an array
@@ -665,6 +660,24 @@ def build_element(self) -> None:
         self.subtype = ", ".join(
             set(component.subtype for component in self.components)
         )
+        if ", " in self.type:
+            warn(
+                f"Variable '{self.name}' is defined with different types:"
+                f" '{self.type}'. This may cause bugs when trying to"
+                "change its value or applying other methods from the "
+                "pysd.py_backend.model.Model class. Running the model "
+                "without modifying this variable should not cause any "
+                "bug."
+            )
+        elif ", " in self.subtype:
+            warn(
+                f"Variable '{self.name}' is defined with different subtypes:"
+                f" '{self.subtype}'. This may cause bugs when trying to"
+                "change its value or applying other methods from the "
+                "pysd.py_backend.model.Model class. Running the model "
+                "without modifying this variable should not cause any "
+                "bug."
+            )
 
     def _manage_multi_def(self, expression: dict) -> str:
         """
@@ -727,13 +740,7 @@ def _build_element_out(self) -> str:
 
         # Get the objects to create as string
         objects = "\n\n".join([
-            value["expression"] % {
-                "final_subs":
-                self.section.subscripts.simplify_subscript_input(
-                    value.get("final_subs", {}))[1]
-            }   # Replace the final subs in the objects that merge
-                # several components
-            for value in self.objects.values()
+            value["expression"] for value in self.objects.values()
             if value["expression"] is not None
         ])
 
diff --git a/pysd/builders/python/subscripts.py b/pysd/builders/python/subscripts.py
index 8bc0c371..62fe8e92 100644
--- a/pysd/builders/python/subscripts.py
+++ b/pysd/builders/python/subscripts.py
@@ -263,7 +263,7 @@ def make_merge_list(self, subs_list: List[List[str]],
                         element: str = "") -> List[str]:
         """
         This is for assisting when building xrmerge. From a list of subscript
-        lists returns the final subscript list after mergin. Necessary when
+        lists returns the final subscript list after merging. Necessary when
         merging variables with subscripts comming from different definitions.
 
         Parameters
diff --git a/pysd/py_backend/external.py b/pysd/py_backend/external.py
index c5930da7..57b2d7bf 100644
--- a/pysd/py_backend/external.py
+++ b/pysd/py_backend/external.py
@@ -748,8 +748,8 @@ def initialize(self):
         """
         Initialize all elements and create the self.data xarray.DataArray
         """
-        if len(self.coordss) == 1:
-            # Just loag one value (no add)
+        if not self.coordss[0]:
+            # Just load one value (no add)
             for self.file, self.sheet, self.x_row_or_col,\
                 self.cell, self.coords\
                 in zip(self.files, self.sheets, self.time_row_or_cols,
@@ -817,7 +817,7 @@ def initialize(self):
         """
         Initialize all elements and create the self.data xarray.DataArray
         """
-        if len(self.coordss) == 1:
+        if not self.coordss[0]:
             # Just loag one value (no add)
             for self.file, self.sheet, self.x_row_or_col,\
                 self.cell, self.coords\
@@ -882,7 +882,7 @@ def initialize(self):
         """
         Initialize all elements and create the self.data xarray.DataArray
         """
-        if len(self.coordss) == 1:
+        if not self.coordss[0]:
             # Just loag one value (no add)
             for self.file, self.sheet, self.transpose, self.cell, self.coords\
                 in zip(self.files, self.sheets, self.transposes,
diff --git a/pysd/py_backend/lookups.py b/pysd/py_backend/lookups.py
index 4f86571a..2624f429 100644
--- a/pysd/py_backend/lookups.py
+++ b/pysd/py_backend/lookups.py
@@ -144,7 +144,18 @@ def initialize(self):
         if len(self.values) == 1:
             # Just loag one value (no add)
             for x, y, coords in self.values:
-                y = np.array(y).reshape((len(x),) + (1,)*len(coords))
+                if len(x) != len(set(x)):
+                    raise ValueError(
+                        self.py_name + "\n"
+                        "x dimension has repeated values..."
+                    )
+                try:
+                    y = np.array(y).reshape((len(x),) + (1,)*len(coords))
+                except ValueError:
+                    raise ValueError(
+                        self.py_name + "\n"
+                        "x and y dimensions have different length..."
+                    )
                 self.data = xr.DataArray(
                     np.tile(y, [1] + utils.compute_shape(coords)),
                     {"lookup_dim": x, **coords},
diff --git a/pysd/py_backend/model.py b/pysd/py_backend/model.py
index d9166cb3..29b36bff 100644
--- a/pysd/py_backend/model.py
+++ b/pysd/py_backend/model.py
@@ -5,6 +5,7 @@
 allows integrating a model or a Macro expression (set of functions in
 a separate file).
 """
+import time
 import warnings
 import inspect
 import pickle
@@ -19,7 +20,8 @@
 
 from . import utils
 from .statefuls import DynamicStateful, Stateful
-from .external import External, Excels
+from .external import External, Excels, ExtLookup, ExtData
+
 from .cache import Cache, constant_cache
 from .data import TabData
 from .lookups import HardcodedLookups
@@ -57,8 +59,8 @@ class Macro(DynamicStateful):
         Time to set at the begginning of the Macro. Default is None.
     data_files: dict or list or str or None
         The dictionary with keys the name of file and variables to
-        load the data from there. Or the list of names or name of the
-        file to search the data in. Only works for TabData type object
+        load the data from. Or the list of names or name of the file
+        to search the data in. Only works for TabData type object
         and it is neccessary to provide it. Default is None.
     py_name: str or None
         The name of the Macro object. Default is None.
@@ -453,13 +455,7 @@ def initialize(self):
 
         if not self.external_loaded:
             # Initialize external elements
-            for element in self._external_elements:
-                element.initialize()
-
-            # Remove Excel data from memory
-            Excels.clean()
-
-            self.external_loaded = True
+            self.initialize_external_data()
 
         # Initialize stateful objects
         for element_name in self.initialize_order:
@@ -532,6 +528,291 @@ def import_pickle(self, file_name):
         self.set_stateful(stateful_dict)
         self.time.set_control_vars(initial_time=time)
 
+    def initialize_external_data(self, externals=None):
+
+        """
+        Initializes external data.
+
+        If a path to a netCDF file containing serialized values of some or
+        all of the model external data is passed in the external argument,
+        those will be loaded from the file.
+
+        To get the full performance gain of loading the externals from a netCDF
+        file, the model should be loaded with initialize=False first.
+
+        Parameters
+        ----------
+        externals: str or pathlib.Path (optional)
+            Path to the netCDF file that contains the model external objects.
+
+        Returns
+        -------
+        None
+
+        """
+
+        if not externals:
+            for ext in self._external_elements:
+                ext.initialize()
+
+            # Remove Excel data from memory
+            Excels.clean()
+
+            self.external_loaded = True
+
+            return
+
+        externals = Path(externals)
+
+        if not externals.is_file():
+            raise FileNotFoundError(f"Invalid file path ({str(externals)})")
+
+        ds = xr.open_dataset(externals)
+
+        for ext in self._external_elements:
+            if ext.py_name in ds.data_vars.keys():
+                # Initialize external from nc file
+                da = ds.data_vars[ext.py_name]
+                if isinstance(ext, ExtData):
+                    # Rename again time dimension
+                    time_dim = [dim for dim in da.dims if dim.startswith(
+                        "time_#")][0]
+                    da = da.rename({time_dim: "time"})
+                elif isinstance(ext, ExtLookup):
+                    # Rename again lookup_dim dimension
+                    lookup_dim = [dim for dim in da.dims if dim.startswith(
+                        "lookup_dim_#")][0]
+                    da = da.rename({lookup_dim: "lookup_dim"})
+                # Assign the value
+                if da.dims:
+                    ext.data = da
+                else:
+                    ext.data = float(da.data)
+            else:
+                # Initialize external from original file
+                ext.initialize()
+
+        Excels.clean()
+
+        self.external_loaded = True
+
+    def serialize_externals(self, export_path="externals.nc",
+                            include_externals="all", exclude_externals=None):
+        """
+        Stores a netCDF file with the data and metadata for all model external
+        objects.
+
+        This method is useful for models with lots of external inputs, which
+        are slow to load into memory. Once exported, the resulting netCDF file
+        can be passed as argument to the initialize_external_data method.
+
+        Names of variables should be those in the model (python safe,
+        without the _ext_type_ string in front).
+
+        Parameters
+        ----------
+        export_path: str or pathlib.Path (optional)
+            Path of the resulting *.nc* file.
+        include_externals: list or str (optional)
+            External objects to export to netCDF.
+            If 'all', then all externals are exported to the *.nc* file.
+            The argument also accepts a list containing spreadsheet file names,
+            external variable names or a combination of both. If a spreadsheet
+            file path is passed, all external objects defined in it will be
+            included in the *.nc* file. Spreadsheet tab names are not currently
+            supported, because the same may be used in different files.
+            Better customisation can be achieved by combining the
+            include_externals and exclude_externals (see description below)
+            arguments.
+        exclude_externals: list or None (optional)
+            Exclude external objects from being included in the exported nc
+            file. It accepts either variable names, spreadsheet files or a
+            combination of both.
+
+        Returns
+        -------
+        None
+
+        """
+        data = {}
+        metadata = {}
+
+        lookup_dims = utils.UniqueDims("lookup_dim")
+        data_dims = utils.UniqueDims("time")
+
+        if isinstance(export_path, str):
+            export_path = Path(export_path)
+
+        if not include_externals:
+            raise ValueError("include_externals argument must not be None.")
+
+        # Generate a Dataframe to make simpler the search of external
+        # objects to include
+        # TODO include also checking the original name
+        py_names = []
+        externals_dict = {"py_var_name": [], "file": [], "ext": []}
+        for ext in self._external_elements:
+            py_names.append(ext.py_name)
+            externals_dict["py_var_name"].append(
+                self.__get_varname_from_ext_name(ext.py_name))
+            externals_dict["file"].append(set(ext.files))
+            externals_dict["ext"].append(ext)
+        exts_df = pd.DataFrame(index=py_names, data=externals_dict)
+
+        if include_externals != "all":
+            if not isinstance(include_externals, (list, set)):
+                raise TypeError(
+                    "include_externals must be 'all', or a list, or a set.")
+
+            # subset only to the externals to include
+            exts_df = exts_df[[
+                name in include_externals
+                or var_name in include_externals
+                or bool(file.intersection(include_externals))
+                for name, (var_name, file)
+                in exts_df[["py_var_name", "file"]].iterrows()
+            ]]
+
+        if exclude_externals:
+            if not isinstance(exclude_externals, (list, set)):
+                raise TypeError("exclude_externals must be a list or a set.")
+
+            # subset only to the externals to include
+            exts_df = exts_df[[
+                name not in exclude_externals
+                and var_name not in exclude_externals
+                and not bool(file.intersection(exclude_externals))
+                for name, (var_name, file)
+                in exts_df[["py_var_name", "file"]].iterrows()
+            ]]
+
+        for _, (ext, var_name) in exts_df[["ext", "py_var_name"]].iterrows():
+            self.__include_for_serialization(
+                ext, var_name, data, metadata, lookup_dims, data_dims
+            )
+
+        # create description to be used as global attribute of the dataset
+        description = {
+            "description": f"External objects for {self.py_model_file} "
+                           f"exported on {time.ctime(time.time())} "
+                           f"using PySD version {__version__}"
+        }
+
+        # create Dataset
+        ds = xr.Dataset(data_vars=data, attrs=description)
+
+        # add data_vars attributes
+        for key, values in metadata.items():
+            ds[key].attrs = values
+
+        ds.to_netcdf(export_path)
+
+    def __include_for_serialization(self, ext, py_name_clean, data, metadata,
+                                    lookup_dims, data_dims):
+        """
+        Initialize the external object and get the data and metadata for
+        inclusion in the netCDF.
+
+        This function updates the metadata dict with the metadata corresponding
+        to each external object, which is collected from model.doc. It also
+        updates the data dict, with the data from the external object (ext).
+
+        It renames the "time" dimension of ExtData types by appending _#
+        followed by a unique number at the end of it. For large models, this
+        prevents having an unnecessary large time dimension, which then causes
+        all ExtData objects to have many nans when stored in a xarray Dataset.
+        It does the same for the "lookup_dim" of all ExtLookup objects.
+
+        NOTE: Though subscripts can be read from Excel, they are hardcoded
+        during the model building process. Therefore they will not be
+        serialized.
+
+        Parameters
+        ----------
+        ext: pysd.py_backend.externals.External
+            External object. It can be any of the External subclasses
+            (ExtConstant, ExtData, ExtLookup)
+        py_name_clean: str
+            Name of the variable without _ext_[constant|data|lookup] prefix.
+        data: dict
+            Collects all the data for each external, which is later used to
+            build the xarray Dataset.
+        metadata: dict
+            Collects the metadata for each external, which is later included
+            as data_vars attributes in the xarray Dataset.
+        lookup_dims: utils.UniqueDims
+            UniqueDims object for "lookup_dim" dimension.
+        data_dims: utils.UniqueDims
+            UniqueDims object for "time" dimension.
+
+        Returns
+        -------
+        None
+
+        """
+        ext.initialize()
+
+        # collecting variable metadata from model._doc
+        var_meta = {
+            col:
+            self._doc.loc[self._doc["Py Name"] == py_name_clean, col].values[0]
+            or "Missing"
+            for col in self._doc.columns
+        }
+        var_meta["files"] = ";".join(ext.files)
+        var_meta["sheets"] = ";".join(ext.sheets)
+        var_meta["cells"] = ";".join(ext.cells)
+        # TODO: add also time_row_or_cols
+
+        da = ext.data
+
+        # Renaming shared dims by all ExtData ("time") and ExtLookup
+        # ("lookup_dim") external objects
+        if isinstance(ext, ExtData):
+            new_name = data_dims.name_new_dim("time", da.coords["time"].values)
+            da = da.rename({"time": new_name})
+        if isinstance(ext, ExtLookup):
+            new_name = lookup_dims.name_new_dim("lookup_dim",
+                                                da.coords["lookup_dim"].values)
+            da = da.rename({"lookup_dim": new_name})
+
+        metadata.update({ext.py_name: var_meta})
+        data.update({ext.py_name: da})
+
+        # TODO use a logger
+        print(f"Finished processing variable {py_name_clean}.")
+
+    def __get_varname_from_ext_name(self, varname):
+        """
+        Returns the name of the variable that depends on the external object
+        named varname. If that is not possible (see warning in the code to
+        understand when that may happen), it gets the name by removing the
+        _ext_[constant|data|lookup] prefix from the varname.
+
+        Parameters
+        ----------
+        varname: str
+            Variable name to which the External object is assigned.
+
+        Returns
+        -------
+        var: str
+            Name of the variable that calls the variable with name varname.
+
+        """
+        for var, deps in self._dependencies.items():
+            for _, ext_name in deps.items():
+                if varname == ext_name:
+                    return var
+
+        warnings.warn(
+            f"No variable depends upon '{varname}'. This is likely due "
+            f"to the fact that '{varname}' is defined using a mix of "
+            "DATA and CONSTANT. Though Vensim allows it, it is "
+            "not recommended."
+        )
+        return "_".join(varname.split("_")[3:])
+
     def get_args(self, param):
         """
         Returns the arguments of a model element.
@@ -1287,18 +1568,18 @@ def select_submodel(self, vars=[], modules=[], exogenous_components={}):
         UserWarning: Selecting submodel, to run the full model again use model.reload()
 
         """
-        c_vars, d_vars, s_deps = self._get_dependencies(vars, modules)
+        deps = self.get_dependencies(vars, modules)
         warnings.warn(
             "Selecting submodel, "
             "to run the full model again use model.reload()")
 
         # get set of all dependencies and all variables to select
-        all_deps = d_vars["initial"].copy()
-        all_deps.update(d_vars["step"])
-        all_deps.update(d_vars["lookup"])
+        all_deps = deps.d_deps["initial"].copy()
+        all_deps.update(deps.d_deps["step"])
+        all_deps.update(deps.d_deps["lookup"])
 
         all_vars = all_deps.copy()
-        all_vars.update(c_vars)
+        all_vars.update(deps.c_vars)
 
         # clean dependendies and namespace dictionaries, and remove
         # the rows from the documentation
@@ -1312,27 +1593,27 @@ def select_submodel(self, vars=[], modules=[], exogenous_components={}):
                 )
 
         for py_name in self._dependencies.copy().keys():
-            if py_name.startswith("_") and py_name not in s_deps:
+            if py_name.startswith("_") and py_name not in deps.s_deps:
                 del self._dependencies[py_name]
 
         # remove active initial from s_deps as they are "fake" objects
         # in dependencies
-        s_deps = {
-            dep for dep in s_deps if not dep.startswith("_active_initial")
+        deps.s_deps = {
+            dep for dep in deps.s_deps if not dep.startswith("_active_initial")
         }
 
         # reassing the dictionary and lists of needed stateful objects
         self._stateful_elements = {
             name: getattr(self.components, name)
-            for name in s_deps
+            for name in deps.s_deps
             if isinstance(getattr(self.components, name), Stateful)
         }
         self._dynamicstateful_elements = [
-            getattr(self.components, name) for name in s_deps
+            getattr(self.components, name) for name in deps.s_deps
             if isinstance(getattr(self.components, name), DynamicStateful)
         ]
         self._macro_elements = [
-            getattr(self.components, name) for name in s_deps
+            getattr(self.components, name) for name in deps.s_deps
             if isinstance(getattr(self.components, name), Macro)
         ]
 
@@ -1392,8 +1673,8 @@ def get_dependencies(self, vars=[], modules=[]):
 
         Returns
         -------
-        dependencies: set
-            Set of dependencies nedded to run vars.
+        dependencies: pysd.py_backend.utils.Dependencies
+            Dependencies data object.
 
         Notes
         -----
@@ -1402,15 +1683,15 @@ def get_dependencies(self, vars=[], modules=[]):
 
         Examples
         --------
-        >>> model.get_dependencies(
-        ...     vars=["Room Temperature", "Teacup temperature"])
+        >>> print(model.get_dependencies(
+        ...     vars=["Room Temperature", "Teacup temperature"]))
         Selected variables (total 1):
             room_temperature, teacup_temperature
         Stateful objects integrated with the selected variables (total 1):
             _integ_teacup_temperature
 
-        >>> model.get_dependencies(
-        ...     modules=["view_1", "view_2/subview_1"])
+        >>> print(model.get_dependencies(
+        ...     modules=["view_1", "view_2/subview_1"]))
         Selected variables (total 4):
             var1, var2, stock1, delay1
         Dependencies for initialization only (total 1):
@@ -1420,9 +1701,9 @@ def get_dependencies(self, vars=[], modules=[]):
         Stateful objects integrated with the selected variables (total 1):
             _integ_stock1, _delay_fixed_delay1
 
-        >>> model.get_dependencies(
+        >>> print(model.get_dependencies(
         ...     vars=["stock3"],
-        ...     modules=["view_1", "view_2/subview_1"])
+        ...     modules=["view_1", "view_2/subview_1"]))
         Selected variables (total 4):
             var1, var2, stock1, stock3, delay1
         Dependencies for initialization only (total 1):
@@ -1431,95 +1712,46 @@ def get_dependencies(self, vars=[], modules=[]):
             _integ_stock1, _integ_stock3, _delay_fixed_delay1
 
         """
-        c_vars, d_vars, s_deps = self._get_dependencies(vars, modules)
-
-        text = utils.print_objects_format(c_vars, "Selected variables")
-
-        if d_vars["initial"]:
-            text += utils.print_objects_format(
-                d_vars["initial"],
-                "\nDependencies for initialization only")
-        if d_vars["step"]:
-            text += utils.print_objects_format(
-                d_vars["step"],
-                "\nDependencies that may change over time")
-        if d_vars["lookup"]:
-            text += utils.print_objects_format(
-                d_vars["lookup"],
-                "\nLookup table dependencies")
-
-        text += utils.print_objects_format(
-            s_deps,
-            "\nStateful objects integrated with the selected variables")
-
-        print(text)
-
-    def _get_dependencies(self, vars=[], modules=[]):
-        """
-        Get the dependencies of a set of variables or modules.
-
-        Parameters
-        ----------
-        vars: set or list of strings (optional)
-            Variables to get the dependencies from.
-            It can be an empty list if the dependencies are computed only
-            using modules. Default is an empty list.
-        modules: set or list of strings (optional)
-            Modules to get the dependencies from.
-            It can be an empty list if the dependencies are computed only
-            using variables. Default is an empty list. Can select a full
-            module or a submodule by passing the path without the .py, e.g.:
-            "view_1/submodule1".
-
-        Returns
-        -------
-        c_vars: set
-            Set of all selected model variables.
-        d_deps: dict of sets
-            Dictionary of dependencies nedded to run vars and modules.
-        s_deps: set
-            Set of stateful objects to update when integrating selected
-            model variables.
-
-        """
-        def check_dep(dependencies, initial=False):
-            for dep in dependencies:
-                if dep in c_vars or dep.startswith("__"):
+        def check_dep(deps_obj, deps, initial=False):
+            for dep in deps:
+                if dep in deps_obj.c_vars or dep.startswith("__"):
                     pass
                 elif dep.startswith("_"):
-                    s_deps.add(dep)
+                    deps_obj.s_deps.add(dep)
                     dep = self._dependencies[dep]
-                    check_dep(dep["initial"], True)
-                    check_dep(dep["step"])
+                    check_dep(deps_obj, dep["initial"], True)
+                    check_dep(deps_obj, dep["step"])
                 else:
-                    if initial and dep not in d_deps["step"]\
-                       and dep not in d_deps["lookup"]:
-                        d_deps["initial"].add(dep)
+                    if initial and dep not in deps_obj.d_deps["step"]\
+                       and dep not in deps_obj.d_deps["lookup"]:
+                        deps_obj.d_deps["initial"].add(dep)
                     else:
-                        if dep in d_deps["initial"]:
-                            d_deps["initial"].remove(dep)
+                        if dep in deps_obj.d_deps["initial"]:
+                            deps_obj.d_deps["initial"].remove(dep)
                         if self.get_args(dep):
-                            d_deps["lookup"].add(dep)
+                            deps_obj.d_deps["lookup"].add(dep)
                         else:
-                            d_deps["step"].add(dep)
+                            deps_obj.d_deps["step"].add(dep)
 
-        d_deps = {"initial": set(), "step": set(), "lookup": set()}
-        s_deps = set()
-        c_vars = {"time", "time_step", "initial_time", "final_time", "saveper"}
+        dependencies = utils.Dependencies(
+            {"time", "time_step", "initial_time", "final_time", "saveper"},
+            {"initial": set(), "step": set(), "lookup": set()},
+            set()
+        )
         for var in vars:
             py_name = utils.get_key_and_value_by_insensitive_key_or_value(
                     var,
                     self._namespace)[1]
-            c_vars.add(py_name)
+            dependencies.c_vars.add(py_name)
         for module in modules:
-            c_vars.update(self.get_vars_in_module(module))
+            dependencies.c_vars.update(self.get_vars_in_module(module))
 
-        for var in c_vars:
+        for var in dependencies.c_vars:
             if var == "time":
                 continue
-            check_dep(self._dependencies[var])
+            check_dep(dependencies, self._dependencies[var])
 
-        return c_vars, d_deps, s_deps
+        return dependencies
 
     def get_vars_in_module(self, module):
         """
diff --git a/pysd/py_backend/output.py b/pysd/py_backend/output.py
index 042e4802..47036896 100644
--- a/pysd/py_backend/output.py
+++ b/pysd/py_backend/output.py
@@ -17,6 +17,7 @@
 import pandas as pd
 
 from pysd._version import __version__
+from pysd.tools.ncfiles import NCFile
 
 from . utils import xrsplit
 
@@ -465,38 +466,10 @@ def postprocess(self, **kwargs):
             self.ds, kwargs["return_addresses"], flatten
             )
         if self.out_file:
-            self.__save_to_file(df)
+            NCFile.df_to_text_file(df, self.out_file)
 
         return df
 
-    def __save_to_file(self, output):
-        """
-        Saves models output.
-
-        Paramters
-        ---------
-        output: pandas.DataFrame
-
-        options: argparse.Namespace
-
-        Returns
-        -------
-        None
-
-        """
-        if self.out_file.suffix == ".tab":
-            sep = "\t"
-        else:
-            sep = ","
-            output.columns = [col.replace(",", ";") for col in output.columns]
-
-        # QUOTE_NONE used to print the csv/tab files as vensim does with
-        # special characterse, e.g.: "my-var"[Dimension]
-        output.to_csv(
-            self.out_file, sep, index_label="Time", quoting=QUOTE_NONE)
-
-        print(f"Data saved in '{self.out_file}'")
-
     def add_run_elements(self, model, capture_elements):
         """
         Adds constant elements to a dataframe.
diff --git a/pysd/py_backend/utils.py b/pysd/py_backend/utils.py
index a66768b7..624f5806 100644
--- a/pysd/py_backend/utils.py
+++ b/pysd/py_backend/utils.py
@@ -8,6 +8,8 @@
 from datetime import datetime
 from pathlib import Path
 from chardet.universaldetector import UniversalDetector
+from dataclasses import dataclass
+from typing import Dict, Set
 
 import progressbar
 import numpy as np
@@ -412,6 +414,49 @@ def print_objects_format(object_set, text):
     return text
 
 
+@dataclass
+class Dependencies():
+    """
+    Representation of variables dependencies.
+
+    Parameters
+    ----------
+    c_vars: set
+        Set of all selected model variables.
+    d_deps: dict
+        Dictionary of dependencies needed to run vars and modules.
+    s_deps: set
+        Set of stateful objects to update when integrating selected
+        model variables.
+
+    """
+    c_vars: Set[str]
+    d_deps: Dict[str, set]
+    s_deps: Set[str]
+
+    def __str__(self):
+        text = print_objects_format(self.c_vars, "Selected variables")
+
+        if self.d_deps["initial"]:
+            text += print_objects_format(
+                self.d_deps["initial"],
+                "\nDependencies for initialization only")
+        if self.d_deps["step"]:
+            text += print_objects_format(
+                self.d_deps["step"],
+                "\nDependencies that may change over time")
+        if self.d_deps["lookup"]:
+            text += print_objects_format(
+                self.d_deps["lookup"],
+                "\nLookup table dependencies")
+
+        text += print_objects_format(
+            self.s_deps,
+            "\nStateful objects integrated with the selected variables")
+
+        return text
+
+
 class ProgressBar:
     """
     Progress bar for integration
@@ -453,3 +498,66 @@ def finish(self):
         except AttributeError:
             # Error if bar is not imported
             pass
+
+
+class UniqueDims():
+    """
+    Helper class to create unique dimension names for data_vars with the
+    same dimension name but different coords in xarray Datasets.
+    """
+    def __init__(self, original_dim_name):
+        self.dim_name = original_dim_name
+        self.dim_prefix = self.dim_name + "_#"
+        self.unique_dims = []
+        self.num = 1
+
+    def name_new_dim(self, dim_name, coords):
+        """
+        Returns either a new name (original_dim_name + _# + num) if the coords
+        list is not in unique_dims, or the preexisting dimension name if it is.
+        Parameters
+        ----------
+        dim_name: str
+            This argument is used to verify that we are passing the right
+            dimension name to the class.
+        coords: list
+            List of coordinates of a dimension.
+
+        Returns
+        -------
+        Updated name of the original dimension.
+        """
+        if dim_name != self.dim_name:
+            raise ValueError(
+                "This object is configured to process dimension "
+                f"{self.dim_name} and it's being passed a dimension "
+                f"named {dim_name}"
+            )
+        if self.is_new(coords):
+            new_dim_name = self.dim_prefix + str(self.num)
+            self.unique_dims.append((new_dim_name, coords))
+            self.num += 1
+            return new_dim_name
+        else:
+            for y in self.unique_dims:
+                if np.array_equal(coords, y[1]):
+                    return y[0]
+
+    def is_new(self, coords):
+        """
+        Checks if coords is already in the unique_dims list or not.
+
+        Parameters
+        ----------
+        coords: list
+            List of coordinates of a dimension.
+
+        Returns
+        -------
+        bool
+        """
+        if not any(
+             map(lambda y: np.array_equal(y[1], coords),
+                 self.unique_dims)):
+            return True
+        return False
diff --git a/pysd/pysd.py b/pysd/pysd.py
index fc032d6b..1c0f9467 100644
--- a/pysd/pysd.py
+++ b/pysd/pysd.py
@@ -42,7 +42,7 @@ def read_xmile(xmile_file, data_files=None, initialize=True,
         If given the list of files where the necessary data to run the model
         is given. Default is None.
 
-    missing_values : str ("warning", "error", "ignore", "keep") (optional)
+    missing_values: str ("warning", "error", "ignore", "keep") (optional)
         What to do with missing values. If "warning" (default)
         shows a warning message and interpolates the values.
         If "raise" raises an error. If "ignore" interpolates
diff --git a/pysd/tools/ncfiles.py b/pysd/tools/ncfiles.py
new file mode 100644
index 00000000..0b2a4939
--- /dev/null
+++ b/pysd/tools/ncfiles.py
@@ -0,0 +1,416 @@
+"""
+Tools for importing and converting netCDF files generated from simulations run
+using PySD.
+"""
+import itertools
+import warnings
+
+from pathlib import Path
+from csv import QUOTE_NONE
+from typing import Union, Optional
+
+import xarray as xr
+import pandas as pd
+
+
+class NCFile():
+    """
+    Helper class to extract data from netCDF files.
+
+    Parameters
+    ----------
+    ncfile: str or pathlib.Path
+        Path to the netCDF file to process.
+    parallel: bool (optional)
+        When True, the Dataset is opened using chunks=-1 (see xarray
+        documentation for details) and DataArrays are processed in parallel
+        using dask delayed. Dask is not included as a requirement for pysd,
+        hence it must be installed separately. Setting parallel=True is
+        highly recommended when the Dataset contains large multidimensional
+        DataArrays.
+
+    """
+
+    valid_export_file_types = [".csv", ".tab"]
+
+    def __init__(self,
+                 filename: Union[str, Path],
+                 parallel: Optional[bool] = False) -> None:
+
+        self.ncfile = NCFile._validate_nc_path(filename)
+        self.parallel = parallel
+        self.ds = self.open_nc()
+
+    def to_text_file(self,
+                     outfile: Optional[Union[str, Path]] = "result.tab",
+                     subset: Optional[list] = None,
+                     time_in_row: Optional[bool] = False,
+                     ) -> pd.DataFrame:
+        """
+        Convert netCDF file contents into comma separated or tab delimited
+        file.
+
+        Parameters
+        ----------
+        outfile: str or pathlib.Path (optional)
+            Path to the output file.
+        subset: list (optional)
+            List of variables to export from the netCDF.
+        time_in_row: bool (optional)
+            Whether time increases along row.
+            Default is False.
+
+        Returns
+        -------
+        df: pandas.DataFrame
+            Dataframe with all colums specified in subset.
+
+        """
+        df = self.to_df(subset=subset)
+
+        NCFile.df_to_text_file(df, outfile, time_in_row)
+
+        return df
+
+    def to_df(self,
+              subset: Optional[list] = None,
+              ) -> pd.DataFrame:
+        """
+        Wrapper to ds_to_df static method. Convert xarray.Dataset into a
+        pandas DataFrame.
+
+        Parameters
+        ----------
+        subset: list (optional)
+            List of variables to export from the Dataset.
+
+        Returns
+        -------
+        df: pandas.DataFrame
+            Dataframe with all colums specified in subset.
+        """
+
+        return NCFile.ds_to_df(self.ds, subset, self.parallel)
+
+    def open_nc(self) -> xr.Dataset:
+        """
+        Loads netCDF file into xarray Dataset. It's basically a wrapper to
+        xr.open_dataset to simplify the interface for pysd use case (loading
+        simulation results).
+
+        Returns
+        -------
+        xarray.Dataset
+
+        """
+        if self.parallel:
+            return xr.open_dataset(self.ncfile, engine="netcdf4", chunks=-1)
+
+        return xr.open_dataset(self.ncfile, engine="netcdf4")
+
+    @staticmethod
+    def ds_to_df(ds: xr.Dataset,
+                 subset: Optional[list] = None,
+                 parallel: Optional[bool] = False,
+                 index_dim: Optional[str] = "time"
+                 ) -> pd.DataFrame:
+        """
+        Convert xarray.Dataset into a pandas DataFrame.
+
+        Parameters
+        ----------
+        ds: xarray.Dataset
+            Dataset object.
+        subset: list (optional)
+            List of variables to export from the Dataset.
+        parallel: bool (optional)
+            When True, DataArrays are processed in parallel using dask delayed.
+            Setting parallel=True is highly recommended when DataArrays are
+            large and multidimensional.
+        index_dim: str (optional)
+            Name of dimensions to use as index of the resulting DataFrame
+            (usually "time").
+
+        Returns
+        -------
+        df: pandas.DataFrame
+            Dataframe with all colums specified in subset.
+
+        """
+        subset = NCFile._validate_ds_subset(ds, subset)
+
+        if parallel:
+            processing_func = NCFile.da_to_dict_delayed
+        else:
+            processing_func = NCFile.da_to_dict
+
+        savedict = {}
+
+        for name in subset:
+            print(f"\nProcessing variable {name}.")
+            da = ds[name]
+            dims = da.dims
+
+            if not dims or dims == (index_dim,):
+                savedict.update({name: da.values.tolist()})
+            else:
+                savedict.update(processing_func(da, index_dim))
+
+        return NCFile.dict_to_df(savedict)
+
+    @staticmethod
+    def df_to_text_file(df: pd.DataFrame, outfile: Path,
+                        time_in_row: Optional[bool] = False
+                        ) -> None:
+        """
+        Store pandas DataFrame into csv or tab file.
+
+        Parameters
+        ----------
+        df: pandas.DataFrame
+            DataFrame to save as csv or tab file.
+        outfile: str or pathlib.Path
+            Path of the output file.
+        time_in_row: bool (optional)
+            Whether time increases along a column or a row.
+
+        Returns
+        -------
+        None
+        """
+        outfile = Path(outfile)
+
+        out_fmt = outfile.suffix
+
+        if out_fmt not in NCFile.valid_export_file_types:
+            raise TypeError("Invalid output file format {out_fmt}\n"
+                            "Supported formats are csv and tab.")
+
+        outfile.parent. mkdir(parents=True, exist_ok=True)
+
+        if not isinstance(time_in_row, bool):
+            raise ValueError("time_in_row argument takes boolen values.")
+
+        # process output file path
+        if outfile.suffix == ".csv":
+            sep = ","
+            df.columns = [col.replace(",", ";") for col in df.columns]
+        else:
+            sep = "\t"
+
+        if time_in_row:
+            df = df.transpose()
+
+        # QUOTE_NONE used to print the csv/tab files as vensim does with
+        # special characterse, e.g.: "my-var"[Dimension]
+        df.to_csv(outfile, sep=sep, index_label="Time", quoting=QUOTE_NONE)
+
+        print(f"Data saved in '{outfile}'")
+
+    @staticmethod
+    def da_to_dict(da: xr.DataArray, index_dim: str) -> dict:
+        """
+        Splits a DataArray into a dictionary, with keys equal to the name
+        of the variable plus all combinations of the cartesian product of
+        coordinates within brackets, and values equal to the data
+        corresponding to those coordinates along the index_dim dimension.
+
+        Parameters
+        ----------
+        index_dim: str
+            The coordinates of this dimension will not be fixed during
+            indexing of the DataArray (i.e. the indexed data will be a
+            scalar or an array along this dimension).
+
+        """
+        dims, coords = NCFile._get_da_dims_coords(da, index_dim)
+
+        indexes = []
+        # TODO: try to achieve the same as itertools.product with
+        # xr.DataArray.stack
+        for coords_prod in itertools.product(*coords):
+            indexes.append(
+                NCFile._index_da_by_coord_labels(da, dims, coords_prod)
+            )
+
+        return dict(indexes)
+
+    @staticmethod
+    def da_to_dict_delayed(da: xr.DataArray, index_dim: str) -> dict:
+        """
+        Same as da_to_dict, but using dask delayed and compute.
+        This function runs much faster when da is a dask array (chunked).
+
+        To use it on its own, you must first make the following imports:
+
+        from dask import delayed, compute
+        from dask.diagnostics import ProgressBar
+
+        Parameters
+        ----------
+        index_dim: str
+            The coordinates of this dimension will not be fixed during
+            indexing (the indexed data will be an array along this dimension).
+
+        """
+        namespace = dir()
+        if not all(
+             map(lambda x: x in namespace, [
+                "delayed", "compute", "ProgressBar"]
+                )):
+            from dask import delayed, compute
+            from dask.diagnostics import ProgressBar
+
+        dims, coords = NCFile._get_da_dims_coords(da, index_dim)
+
+        # loading data into memory for faster indexing
+        da.load()
+
+        indexes = []
+        # TODO: try to achieve the same as itertools.product with
+        # xr.DataArray.stack
+        for coords_prod in itertools.product(*coords):
+            x = delayed(
+                NCFile._index_da_by_coord_labels
+            )(da, dims, coords_prod)
+            indexes.append(x)
+
+        with ProgressBar():
+            res = compute(*indexes)
+
+        return dict(res)
+
+    @staticmethod
+    def dict_to_df(d: dict) -> pd.DataFrame:
+        """
+        Convert a dict to a pandas Dataframe.
+
+        Parameters
+        ----------
+        d: dict
+            Dictionary to convert to pandas DataFrame.
+
+        """
+        if "time" not in d:
+            raise KeyError("Missing time key.")
+
+        return pd.DataFrame(d).set_index('time')
+
+    @staticmethod
+    def _validate_nc_path(nc_path: Union[str, Path]) -> Path:
+        """
+        Checks validity of the nc_path passed by the user. We run these
+        checks because xarray Exceptions are not very explicit.
+        """
+
+        if not isinstance(nc_path, (str, Path)):
+            raise TypeError(f"Invalid file path type: {type(nc_path)}.\n"
+                            "Please provide string or pathlib Path")
+
+        nc_path = Path(nc_path)
+
+        if not nc_path.is_file():
+            raise FileNotFoundError(f"{nc_path} could not be found.")
+
+        if not nc_path.suffix == ".nc":
+            raise ValueError("Input file must have nc extension.")
+
+        return nc_path
+
+    @staticmethod
+    def _validate_ds_subset(ds: xr.Dataset, subset: list) -> list:
+        """
+        If subset=None, it returns a list with all variable names in the ds.
+        If var names in subset are present in ds, it returns them, else it
+        warns the user.
+
+        Parameters
+        ----------
+        subset: list
+            Subset of variable names in the xarray Dataset.
+
+        """
+        # use all variable names
+        if not subset:
+            new_subset = [name for name in ds.data_vars.keys()]
+        else:
+            if not isinstance(subset, list) or \
+                 not all(map(lambda x: isinstance(x, str), subset)):
+                raise TypeError("Subset argument must be a list of strings.")
+
+            new_subset = []
+            for name in subset:
+                if name in ds.data_vars.keys():
+                    new_subset.append(name)
+                else:
+                    warnings.warn(f"{name} not in Dataset.")
+
+            if not new_subset:
+                raise ValueError("None of the elements of the subset are "
+                                 "present in the Dataset.")
+
+        # adding time in the final subset
+        if "time" not in new_subset:
+            new_subset.append("time")
+
+        return new_subset
+
+    @staticmethod
+    def _index_da_by_coord_labels(da: xr.DataArray, dims: list,
+                                  coords: tuple) -> tuple:
+        """
+        Generates variable names, combining the actual name of the variable
+        with the coordinate names between brackets and separated by commas,
+        and indexes the DataArray by the coordinate names specified in
+        the coords argument.
+
+        Parameters
+        ----------
+        da: xr.Dataset
+            Dataset to be indexed.
+        dims: list
+            Dimensions along which the DataArray will be indexed.
+        coords: tuple
+            Coordinate names for each of the dimensons in the dims list.
+
+        Returns
+        -------
+        A tuple consisting of the string
+        var_name[dim_1_coord_j, ..., dim_n_coord_k] in the first index, and the
+        indexed data as the second index.
+
+        """
+        name = da.name
+        idx = dict(zip(dims, coords))
+        subs = "[" + ",".join(map(lambda x: str(x), coords)) + "]"
+
+        return name + subs, da.loc[idx].values
+
+    @staticmethod
+    def _get_da_dims_coords(da: xr.DataArray, exclude_dim: str) -> tuple:
+        """
+        Returns the dimension names and coordinate labels in two
+        separate lists. If a dimension name is in the exclude_dims
+        list, the returned dims and coords will not include it.
+
+        Parameters
+        ----------
+        exclude_dim: str
+            Names of dimension to exclude.
+
+        Returns
+        -------
+        dims: list
+            List containing the names of the DataArray dimensions.
+        coords: list
+            List of lists of coordinates for each dimension.
+
+        """
+        dims, coords = [], []
+
+        for dim in da.dims:
+            if dim != exclude_dim:
+                dims.append(dim)
+                coords.append(da.coords[dim].values)
+
+        return dims, coords


Test patch:
diff --git a/tests/conftest.py b/tests/conftest.py
index 0019e9ad..fe3f39d2 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -22,9 +22,9 @@ def _test_models(_root):
 
 
 @pytest.fixture(scope="class")
-def shared_tmpdir(tmpdir_factory):
+def shared_tmpdir(tmp_path_factory):
     # shared temporary directory for each class
-    return Path(tmpdir_factory.mktemp("shared"))
+    return tmp_path_factory.mktemp("shared")
 
 
 @pytest.fixture
diff --git a/tests/pytest.ini b/tests/pytest.ini
index ac62136e..715bcb12 100644
--- a/tests/pytest.ini
+++ b/tests/pytest.ini
@@ -4,3 +4,6 @@ filterwarnings =
     error
     ignore:Creating an ndarray from ragged nested sequences
     ignore:`np.bool` is a deprecated alias for the builtin `bool`. To silence this warning, use `bool` by itself. Doing this will not modify any behavior and is safe. If you specifically wanted the numpy scalar type, use `np.bool_` here.
+    ignore:SelectableGroups dict interface is deprecated. Use select.:DeprecationWarning
+    ignore:tostring\(\) is deprecated. Use tobytes\(\) instead.:DeprecationWarning
+    ignore:distutils Version classes are deprecated. Use packaging.version instead.:DeprecationWarning
diff --git a/tests/pytest_pysd/pytest_select_submodel.py b/tests/pytest_pysd/pytest_select_submodel.py
index c4759aea..d768db14 100644
--- a/tests/pytest_pysd/pytest_select_submodel.py
+++ b/tests/pytest_pysd/pytest_select_submodel.py
@@ -87,38 +87,38 @@ def test__get_dependencies(self, model, variables, modules,
                                n_deps, dep_vars):
 
         # get selected vars and dependencies in sets and dictionary
-        out = model._get_dependencies(vars=variables, modules=modules)
-        assert len(out[0]) == n_deps[0]
-        assert len(out[1]["initial"]) == n_deps[1]
-        assert len(out[1]["step"]) == n_deps[2]
-        assert len(out[1]["lookup"]) == n_deps[3]
-        assert len(out[2]) == n_deps[4]
+        out = model.get_dependencies(vars=variables, modules=modules)
+        assert len(out.c_vars) == n_deps[0]
+        assert len(out.d_deps["initial"]) == n_deps[1]
+        assert len(out.d_deps["step"]) == n_deps[2]
+        assert len(out.d_deps["lookup"]) == n_deps[3]
+        assert len(out.s_deps) == n_deps[4]
 
-        assert self.common_vars.issubset(out[0])
+        assert self.common_vars.issubset(out.c_vars)
 
     def test_get_dependencies(self, capsys, model, variables, modules,
                               n_deps, dep_vars):
         # get the dependencies information of the selected variables and
-        # modules by stdout
-        model.get_dependencies(vars=variables, modules=modules)
+        # modules as string
+        deps = model.get_dependencies(vars=variables, modules=modules)
 
-        captured = capsys.readouterr()  # capture stdout
+        print_message = deps.__str__()
 
         for n, message in zip(n_deps, self.messages):
             if n != 0:
                 # check the message with the number of dependencies of
                 # each type
-                assert message + " (total %s):\n" % n in captured.out
+                assert message + " (total %s):\n" % n in print_message
             else:
                 # if not dependencies not message should be printed
-                assert message not in captured.out
+                assert message not in print_message
 
         # assert _integ_stock is in the message as the included stateful object
-        assert "_integ_stock" in captured.out
+        assert "_integ_stock" in print_message
 
         # assert all dependencies of the submodel are in the message
         for var in dep_vars:
-            assert var in captured.out
+            assert var in print_message
 
     def test_select_submodel(self, model, variables, modules,
                              n_deps, dep_vars):
diff --git a/tests/pytest_pysd/pytest_serialize_externals.py b/tests/pytest_pysd/pytest_serialize_externals.py
new file mode 100644
index 00000000..33f45d77
--- /dev/null
+++ b/tests/pytest_pysd/pytest_serialize_externals.py
@@ -0,0 +1,351 @@
+import pytest
+
+from pathlib import Path
+import xarray as xr
+import numpy as np
+
+import pysd
+
+data_ext = Path(
+    "test-models/tests/get_data_args_3d_xls/test_get_data_args_3d_xls.mdl")
+constant_ext = Path("test-models/tests/get_constants_subranges/"
+                    "test_get_constants_subranges.mdl")
+lookup_ext = Path("test-models/tests/get_lookups_data_3d_xls/"
+                  "test_get_lookups_data_3d_xls.mdl")
+lookup_subranges = Path("test-models/tests/get_lookups_subset/"
+                        "test_get_lookups_subset.mdl")
+mixed_definitions = Path("test-models/tests/get_mixed_definitions/"
+                         "test_get_mixed_definitions.mdl")
+subscript_ext = Path("test-models/tests/get_subscript_3d_arrays_xls/"
+                     "test_get_subscript_3d_arrays_xls.mdl")
+multiple_excels = Path("test-models/tests/get_lookups_data_3d_xls/"
+                       "test_get_lookups_data_3d_xls.mdl")
+incomplete_constant_def = Path("test-models/tests/"
+                               "get_constants_incomplete_subscript/"
+                               "test_get_constants_incomplete_subscript.mdl")
+
+
+class TestSerialization():
+
+    @pytest.mark.parametrize("model_path", [data_ext])
+    def test_serialize_externals_file_content(self, tmp_path, model):
+
+        result_path = tmp_path / "externals.nc"
+
+        model.serialize_externals(export_path=result_path,
+                                  include_externals=["input.xls"],
+                                  exclude_externals=["data_backward"])
+
+        ds = xr.open_dataset(result_path)
+
+        assert "_ext_data_data_backward" not in ds.data_vars.keys()
+        assert "_ext_data_data_forward" in ds.data_vars.keys()
+        assert np.allclose(
+            ds.data_vars["_ext_data_data_forward"].loc[
+                {"time_#1": 10.0}].data,
+            np.array([[10., 5.], [-5., 1.]]), equal_nan=True)
+        assert np.allclose(ds.coords["time_#1"].values,
+                           np.array([0.,  5., 10., 15., 20., 25., 30.]),
+                           equal_nan=True)
+        assert "description" in ds.attrs.keys()
+        assert all(map(
+            lambda x: x in ds.data_vars["_ext_data_data_forward"].attrs.keys(),
+            ["Py Name", "sheets", "files"]))
+
+    @pytest.mark.parametrize(
+        "model_path,include,exclude",
+        [
+            (data_ext, ["data_forward"], None),
+            (data_ext, "all", ["data_backward"]),
+            (lookup_ext, ["lookup_function"], None),
+            (lookup_subranges, "all", None),
+            (constant_ext, ["constant"], None),
+        ],
+        ids=["data", "data_exclude", "lookup", "lookup_subranges", "constant"])
+    @pytest.mark.filterwarnings("ignore")
+    def test_serialize_externals_different_types(self, model,
+                                                 include, exclude):
+
+        # external subscripts are not included because they are hardcoded in
+        # the Python model.
+        model_path = Path(model.py_model_file)
+        externals_path = model_path.with_suffix(".nc")
+
+        if include == "all":
+            if exclude:
+                varnames_included = {
+                    ext.py_name: "_".join(ext.py_name.split("_")[3:])
+                    for ext in model._external_elements
+                    if not "_".join(ext.py_name.split("_")[3:]) in exclude
+                }
+            else:
+                varnames_included = {
+                    ext.py_name: "_".join(ext.py_name.split("_")[3:])
+                    for ext in model._external_elements
+                }
+        else:
+            if exclude:
+                varnames_included = {
+                    model._dependencies[var]["__external__"]: var
+                    for var in include if var not in exclude
+                }
+            else:
+                varnames_included = {
+                    model._dependencies[var]["__external__"]: var
+                    for var in include
+                }
+
+        model.serialize_externals(export_path=externals_path,
+                                  include_externals=include,
+                                  exclude_externals=exclude)
+
+        ds = xr.open_dataset(externals_path)
+
+        # all expected variables are in the dataset
+        assert all([x in ds.data_vars.keys() for x in varnames_included])
+
+        if exclude:
+            assert not any([x in ds.data_vars.keys() for x in exclude])
+
+        # the values in the dataset are the same than in the model
+        assert all(np.allclose(
+            ds.data_vars[py_name].values,
+            getattr(model.components, py_name).data.values, equal_nan=True)
+                for py_name in varnames_included.keys())
+
+        model = pysd.load(model_path, initialize=False)
+        model.initialize_external_data(externals=externals_path)
+
+        # the values in the dataset are the same than in the model
+        # reinitialized
+        assert all(np.allclose(
+            ds.data_vars[py_name].values,
+            getattr(model.components, py_name).data.values, equal_nan=True)
+                for py_name in varnames_included.keys())
+
+    @pytest.mark.parametrize("model_path", [mixed_definitions])
+    @pytest.mark.filterwarnings("ignore")
+    def test_serialize_mixed_definitions(self, model):
+
+        model_path = Path(model.py_model_file)
+        externals_path = model_path.with_suffix(".nc")
+
+        varname = "_ext_constant_variable"
+        warn_message =\
+            f"No variable depends upon '{varname}'. This is likely due to"\
+            f" the fact that '{varname}' is defined using a mix of DATA"\
+            " and CONSTANT. Though Vensim allows it, it is not recommended."
+        with pytest.warns(UserWarning, match=warn_message):
+            model.serialize_externals(export_path=externals_path,
+                                      include_externals="all",
+                                      exclude_externals=None)
+
+        ds = xr.open_dataset(externals_path)
+
+        assert list(ds.coords) == [
+            'dim', 'adim', 'Type', 'DimAB', 'time_#1', 'time_#2']
+        assert ds.data_vars["_ext_constant_const_var_1"].dims == (
+            'dim',)
+
+        # There will be nans stored in the dataset if the dimension
+        # does not have values in all dimensions (incomplete definition)
+        assert any(ds.data_vars["_ext_constant_const_var_1"].isnull())
+        np.testing.assert_array_equal(
+            model["_ext_constant_const_var_1"].data,
+            ds.data_vars["_ext_constant_const_var_1"].data)
+
+        # re stands for reinitialized loading from external data file
+        model_re = pysd.load(model_path, initialize=False)
+        model_re.initialize_external_data(externals=externals_path)
+        assert model_re.external_loaded
+
+        # Once the model is re initialized loading external data, it should not
+        # have the nans that are in the dataset, even if it is an incompletely
+        # defined index
+        np.testing.assert_array_equal(
+            model_re["_ext_constant_const_var_1"].data,
+            model["_ext_constant_const_var_1"].data)
+        np.testing.assert_array_equal(
+            model_re["_ext_constant_const_var_2"].data,
+            model["_ext_constant_const_var_2"].data)
+        # this one has nans when initializing from excel
+        np.testing.assert_array_equal(
+            model_re["_ext_constant_const_var_3"].data,
+            model["_ext_constant_const_var_3"].data)
+        np.testing.assert_array_equal(
+            model_re["_ext_constant_variable"].data,
+            model["_ext_constant_variable"].data)
+        np.testing.assert_array_equal(
+            model_re.get_series_data("_ext_data_data_var_1").data,
+            model.get_series_data("_ext_data_data_var_1").data)
+        np.testing.assert_array_equal(
+            model_re.get_series_data("_ext_data_data_var_2").data,
+            model.get_series_data("_ext_data_data_var_2").data)
+        np.testing.assert_array_equal(
+            model_re.get_series_data("_ext_data_data_var_3").data,
+            model.get_series_data("_ext_data_data_var_3").data)
+        np.testing.assert_array_equal(
+            model_re.get_series_data("_ext_data_variable").data,
+            model.get_series_data("_ext_data_variable").data)
+        np.testing.assert_array_equal(
+            model_re.get_series_data("_ext_data_variable").data,
+            model.get_series_data("_ext_data_variable").data)
+        # check that the time coordinates are the same in the ds than in the
+        # model
+        np.testing.assert_array_equal(
+            ds.data_vars["_ext_data_variable"].coords["time_#2"],
+            model_re.get_series_data("_ext_data_variable").coords["time"]
+            )
+
+        # data interpolation
+        np.testing.assert_array_equal(
+            getattr(model_re.components, "_ext_data_data_var_3")(1.5).data,
+            getattr(model.components, "_ext_data_data_var_3")(1.5).data
+            )
+        np.testing.assert_array_equal(
+            getattr(model_re.components, "_ext_data_variable")(5.5).data,
+            getattr(model.components, "_ext_data_variable")(5.5).data)
+
+        warn_message = r"extrapolating data (above|below) the "\
+            r"(maximum|minimum) value of the time"
+        with pytest.warns(UserWarning, match=warn_message):
+            np.testing.assert_array_equal(
+                getattr(model_re.components, "_ext_data_variable")(25).data,
+                getattr(model.components, "_ext_data_variable")(25).data)
+            np.testing.assert_array_equal(
+                getattr(model_re.components, "_ext_data_variable")(-1).data,
+                getattr(model.components, "_ext_data_variable")(-1).data)
+
+        # they should also have the same dimension names
+        assert model_re["_ext_constant_const_var_1"].dims == \
+            model["_ext_constant_const_var_1"].dims
+        assert model_re["_ext_constant_const_var_2"].dims == \
+            model["_ext_constant_const_var_2"].dims
+        assert model_re["_ext_constant_const_var_3"].dims == \
+            model["_ext_constant_const_var_3"].dims
+        assert model_re["_ext_constant_variable"].dims == \
+            model["_ext_constant_variable"].dims
+        assert model_re.get_series_data("_ext_data_data_var_1").dims == \
+            model.get_series_data("_ext_data_data_var_1").dims
+        assert model_re.get_series_data("_ext_data_data_var_2").dims == \
+            model.get_series_data("_ext_data_data_var_2").dims
+        assert model_re.get_series_data("_ext_data_data_var_3").dims == \
+            model.get_series_data("_ext_data_data_var_3").dims
+        assert model_re.get_series_data("_ext_data_variable").dims == \
+            model.get_series_data("_ext_data_variable").dims
+
+    @pytest.mark.parametrize("model_path", [subscript_ext])
+    def test_serialize_ext_subscript(self, model, tmp_path):
+
+        # subscripts read from external file sare hardcoded in the python
+        # translation, hence they should not be present in the nc file.
+        result_path = tmp_path / "externals_subs.nc"
+
+        model.serialize_externals(export_path=result_path,
+                                  include_externals="all",
+                                  exclude_externals=None)
+
+        ds = xr.open_dataset(result_path)
+
+        assert len(ds.data_vars.keys()) == 1
+
+    @pytest.mark.parametrize(
+        "model_path,include,exclude,included,excluded",
+        [
+            (
+                multiple_excels, "all", ["input2.xls"],
+                ["_ext_lookup_lookup_function"],
+                ["_ext_data_data_function_table"]
+            ),
+            (
+                multiple_excels, "all", ["input.xls"],
+                ["_ext_data_data_function_table"],
+                ["_ext_lookup_lookup_function"]
+            ),
+            (
+                multiple_excels, "all", None,
+                ["_ext_data_data_function_table",
+                 "_ext_lookup_lookup_function"],
+                []
+            ),
+            (
+                multiple_excels, ["input.xls", ], None,
+                ["_ext_lookup_lookup_function"],
+                []
+            ),
+            (
+                multiple_excels, ["input.xls", "input2.xls"], None,
+                ["_ext_data_data_function_table",
+                 "_ext_lookup_lookup_function"],
+                []
+            ),
+            (
+                multiple_excels, ["input.xls", "data_function_table"], None,
+                ["_ext_data_data_function_table",
+                 "_ext_lookup_lookup_function"],
+                []
+            ),
+         ]
+    )
+    def test_serialize_combine_vars_and_excels(self, model, include, exclude,
+                                               included, excluded):
+
+        # This test combines variable names an excel sheets.
+        model_path = Path(model.py_model_file)
+        externals_path = model_path.with_suffix(".nc")
+
+        model.serialize_externals(export_path=str(externals_path),
+                                  include_externals=include,
+                                  exclude_externals=exclude)
+
+        ds = xr.open_dataset(externals_path)
+
+        assert all(x in ds.data_vars.keys() for x in included)
+        assert not any(x in ds.data_vars.keys() for x in excluded)
+
+    @pytest.mark.parametrize("model_path", [incomplete_constant_def])
+    @pytest.mark.filterwarnings("ignore")
+    def test_incomplete_constant_definition(self, model):
+
+        py_model = Path(model.py_model_file)
+        ext_path = py_model.parent / "externals.nc"
+
+        model.serialize_externals(ext_path, include_externals="all")
+
+        model_re = pysd.load(py_model, initialize=False)
+        model_re.initialize_external_data(ext_path)
+
+        assert model["adimensional_const"] == model_re["adimensional_const"]
+
+        # the reinitialized model will have nans, because the final_coords is
+        # not consistent with the actual coords of the variable
+        np.testing.assert_array_equal(
+            model["_ext_constant_var_only_women"].data,
+            model_re["_ext_constant_var_only_women"].data)
+
+        assert getattr(
+            model.components,
+            "_ext_constant_var_only_women").final_coords["dim2"]\
+            == ["female", "male"]
+
+    @pytest.mark.parametrize("model_path", [subscript_ext])
+    def test_exceptions(self, model):
+
+        # not including any externals should fail
+        with pytest.raises(ValueError):
+            model.serialize_externals(export_path="externals.nc",
+                                      include_externals=None)
+
+        # include externals must be "all" or list
+        with pytest.raises(TypeError):
+            model.serialize_externals(export_path="externals.nc",
+                                      include_externals="hello_world")
+
+        # exclude externals must be a list
+        with pytest.raises(TypeError):
+            model.serialize_externals(export_path="externals.nc",
+                                      include_externals="all",
+                                      exclude_externals="hello_world")
+        # externals file not found
+        with pytest.raises(FileNotFoundError):
+            model.initialize_external_data(externals="missing_file.nc")
diff --git a/tests/pytest_pysd/pytest_utils.py b/tests/pytest_pysd/pytest_utils.py
index a3210c07..221c5771 100644
--- a/tests/pytest_pysd/pytest_utils.py
+++ b/tests/pytest_pysd/pytest_utils.py
@@ -7,7 +7,7 @@
 from pysd.tools.benchmarking import assert_frames_close
 import pysd
 from pysd.py_backend.utils import compute_shape, rearrange, load_outputs,\
-    ProgressBar
+    ProgressBar, UniqueDims
 
 
 class TestUtils():
@@ -217,6 +217,32 @@ def test_rearrange(self):
 
         assert rearrange(None, ['d2'], _subscript_dict) is None
 
+    @pytest.mark.parametrize("dim_name,inputs,expected_unique,expected_return",
+    [
+        ("dim1", [[1, 2, 3], [1, 2, 3]], [("dim1_#1", [1, 2, 3])],
+        ["dim1_#1", "dim1_#1"]),
+        ("dim2", [[1, 2, 3], [2, 3], [3], [1, 2, 3, 4], [1, 2, 3]],
+        [("dim2_#1", [1, 2, 3]), ("dim2_#2", [2, 3]), ("dim2_#3", [3]),
+        ("dim2_#4", [1, 2, 3, 4])],
+        ["dim2_#1", "dim2_#2", "dim2_#3", "dim2_#4", "dim2_#1"]),
+    ]
+    )
+    def test_unique_dims_in_dataset(self, dim_name, inputs, expected_unique,
+                                    expected_return):
+
+        unique_dims = UniqueDims(dim_name)
+
+        for coords, ret in zip(inputs, expected_return):
+            assert unique_dims.name_new_dim(dim_name, coords) == ret
+
+        error_msg = "This object is configured to process dimension "\
+                  + f"{dim_name} and it's being passed a dimension named "\
+                  + "wrong_dim_name"
+
+        with pytest.raises(ValueError, match=error_msg):
+            unique_dims.name_new_dim("wrong_dim_name", coords)
+
+        assert unique_dims.unique_dims == expected_unique
 
 class TestLoadOutputs():
     def test_non_valid_outputs(self, _root):
diff --git a/tests/pytest_tools/pytest_ncfile.py b/tests/pytest_tools/pytest_ncfile.py
new file mode 100644
index 00000000..76c875a2
--- /dev/null
+++ b/tests/pytest_tools/pytest_ncfile.py
@@ -0,0 +1,212 @@
+import pytest
+
+from pathlib import Path
+
+import xarray as xr
+import numpy as np
+import pandas as pd
+
+from pysd.tools.ncfiles import NCFile
+
+@pytest.fixture(scope="session")
+def sample_dataset():
+    # generation of data (slight modification from xarray.Dataset example)
+    np.random.seed(0)
+    temperature = 15 + 8 * np.random.randn(3, 2, 2)
+    precipitation = 10 * np.random.rand(3, 2, 2)
+    altitude = 1000 * np.random.rand(2, 2)
+    lon = [[-99.83, -99.32], [-99.79, -99.23]]
+    lat = [[42.25, 42.21], [42.63, 42.59]]
+    time = pd.date_range("2014-09-06", periods=3)
+
+    ds = xr.Dataset(
+        data_vars=dict(
+            temperature=(
+                ["time", "x", "y"], temperature),
+            precipitation=(
+                ["time", "x", "y"], precipitation),
+            altitude=(["x", "y"], altitude),
+            non_dimensional_var=([], np.array(5)),
+            time=(["time"], time)
+        ),
+        coords=dict(
+            lon=(["x", "y"], lon),
+            lat=(["x", "y"], lat),
+        ),
+        attrs=dict(description="Weather related data."),
+    )
+
+    return ds
+
+@pytest.fixture(scope="session")
+def nc_file(sample_dataset, tmp_path_factory):
+
+    path = tmp_path_factory.mktemp("data") / "dataset.nc"
+
+    sample_dataset.to_netcdf(
+        path, mode="w", format="NETCDF4", engine="netcdf4")
+
+    return path
+
+
+class TestNCFile():
+
+    @pytest.mark.parametrize(
+        "file_path,raises,expect",
+        [(1, True, TypeError),
+        ("unexisting_file.nc", True, FileNotFoundError),
+        (Path(__file__), True, ValueError),
+        ("nc_file", False, "nc_file")
+        ],
+        ids=["Wrong type", "Not found", "Wrong extension", "Correct path"]
+        )
+    def test__validate_nc_path_errors(
+         self, request, file_path, raises, expect):
+        if not raises:
+            nc_file = request.getfixturevalue(file_path)
+            assert NCFile._validate_nc_path(nc_file) == nc_file
+        else:
+            with pytest.raises(expect):
+                NCFile._validate_nc_path(file_path)
+
+    @pytest.mark.parametrize(
+        "subset,raises,expected",
+        [
+            ([], False, [
+                            "temperature",
+                            "precipitation",
+                            "altitude",
+                            "non_dimensional_var",
+                            "time"
+                        ]
+            ),
+            (["temperature"], False, ["temperature", "time"]),
+            (["non_existing"], True, ValueError),
+            ("I am not a list", True, TypeError)
+        ],
+        ids=["empty", "one var", "non existing single var", "wrong type"]
+        )
+    @pytest.mark.filterwarnings("ignore")
+    def test__validate_ds_subset(
+         self, sample_dataset, subset, raises, expected):
+
+        if raises:
+            with pytest.raises(expected):
+                NCFile._validate_ds_subset(sample_dataset, subset)
+        else:
+            assert NCFile._validate_ds_subset(
+                sample_dataset, subset) == expected
+
+    @pytest.mark.parametrize(
+        "data_var,dims,coords,expected",
+        [
+            ("temperature", ["x", "y"], (0,0), ("temperature[0,0]", np.array(
+                [29.11241876774131 ,29.94046392119974, 14.174249185651536]))
+            ),
+            ("altitude", ["x", "y"], (1,1), ("altitude[1,1]", np.array(
+                944.66891705))
+            ),
+        ],
+        ids = ["time-dependent data var", # returns vector
+                "constant data var"] # returns scalar
+        )
+    def test__index_da_by_coord_labels(
+        self, sample_dataset, data_var, dims, coords, expected):
+
+        da = sample_dataset.data_vars[data_var]
+
+        result = NCFile._index_da_by_coord_labels(
+            da, dims, coords)
+
+        assert result[0] == expected[0]
+        assert np.allclose(result[1], expected[1])
+
+    @pytest.mark.parametrize(
+        "data_var,expected_keys",
+        [
+        ("temperature",
+        [
+            'temperature[0,0]',
+            'temperature[0,1]',
+            'temperature[1,0]',
+            'temperature[1,1]'
+        ]),
+        ("altitude",
+        [
+            'altitude[0,0]',
+            'altitude[0,1]',
+            'altitude[1,0]',
+            'altitude[1,1]'
+        ]),
+        ],
+        ids = ["time-dependent data var", "constant data var"]
+        )
+    def test_da_to_dict(self, sample_dataset, data_var, expected_keys):
+        idx = "time"
+        da = sample_dataset.data_vars[data_var]
+
+        serial = NCFile.da_to_dict(da, index_dim=idx)
+        assert all(map(lambda x: x in serial.keys(), expected_keys))
+
+        delayed = NCFile.da_to_dict_delayed(da, index_dim=idx)
+
+        assert all(
+            map(lambda key: np.allclose(serial[key], delayed[key]),
+                expected_keys
+            )
+        )
+
+    @pytest.mark.parametrize(
+        "d,raises,expected",
+        [
+        ({"time": [1, 2, 3], "col1": 4, "col2": [2, 3, 4]}, False,
+        [[1, 2, 3], [4]*3]
+        ),
+        ({"col1": 4, "col2": [2, 3, 4]}, True, "irrelevant")
+        ],
+        ids=["ok", "missing time"])
+    def test_dict_to_df(self, d, raises, expected):
+
+        if raises:
+            with pytest.raises(KeyError):
+                NCFile.dict_to_df(d)
+        else:
+            df = NCFile.dict_to_df(d)
+
+            assert all(df.index == expected[0])
+            assert all(df["col1"] == expected[1])
+
+    def test_df_to_text_file_errors(self, nc_file):
+
+        with pytest.raises(TypeError):
+            NCFile.df_to_text_file(None, str(nc_file), time_in_row=False)
+
+        csv_path = nc_file.parent / (nc_file.name + '.csv')
+        with pytest.raises(ValueError):
+            NCFile.df_to_text_file(None, csv_path, time_in_row="False")
+
+    @pytest.mark.parametrize(
+        "outfmt,time_in_row,parallel", [
+            (".csv", True, True),
+            (".csv", True, False),
+            (".csv", False, True),
+            (".csv", False, False),
+            (".tab", True, True),
+            (".tab", True, False),
+            (".tab", False, True),
+            (".tab", False, False)
+            ])
+    def test_to_text_file(
+         self, shared_tmpdir, nc_file, outfmt, time_in_row, parallel):
+        obj = NCFile(nc_file, parallel=parallel)
+        assert obj.ncfile == nc_file
+        assert isinstance(obj.ds, xr.Dataset)
+
+        outfile = shared_tmpdir / f"data_from_nc{outfmt}"
+        obj.to_text_file(outfile, time_in_row=time_in_row)
+        assert outfile.is_file()
+
+    def test_ds_to_df(self, sample_dataset):
+        df = NCFile.ds_to_df(sample_dataset, subset=["temperature"])
+        assert all(map(lambda x: x.startswith("temperature"),df.columns))
+        assert df.shape == (3, 4)
diff --git a/tests/pytest_types/lookup/pytest_lookups.py b/tests/pytest_types/lookup/pytest_lookups.py
index 490fc025..740c0607 100644
--- a/tests/pytest_types/lookup/pytest_lookups.py
+++ b/tests/pytest_types/lookup/pytest_lookups.py
@@ -3,7 +3,7 @@
 import xarray as xr
 import pandas as pd
 
-from pysd.py_backend.lookups import Lookups
+from pysd.py_backend.lookups import Lookups, HardcodedLookups
 
 
 @pytest.mark.parametrize(
@@ -103,3 +103,38 @@ def test_lookups_set_value(self, lookups, new_value, expected):
             assert lookups.data == expected
         else:
             assert lookups.data.equals(expected)
+
+
+@pytest.mark.parametrize(
+    "x,y,error_type,error_message",
+    [
+        (  # repeated
+            [1, 2, 2, 3, 4],
+            [10, 20, 20, 30, 40],
+            ValueError,
+            r"my_look\nx dimension has repeated values..."
+        ),
+        (  # diff-len1
+            [1, 2, 3, 4],
+            [10, 20, 20, 30, 40],
+            ValueError,
+            r"my_look\nx and y dimensions have different length..."
+        ),
+        (  # diff-len2
+            [1, 2, 3, 4],
+            [10, 30, 40],
+            ValueError,
+            r"my_look\nx and y dimensions have different length..."
+        ),
+    ],
+    ids=[
+        "repeated", "diff-len1", "diff-len2"
+    ]
+)
+class TestLookupsError():
+
+    def test_hardcoded_lookups_error(self, x, y, error_type, error_message):
+        hl = HardcodedLookups(x, y, {}, None, {}, "my_look")
+
+        with pytest.raises(error_type, match=error_message):
+            hl.initialize()
diff --git a/tests/requirements.txt b/tests/requirements.txt
index 833d24f2..6d8ec202 100644
--- a/tests/requirements.txt
+++ b/tests/requirements.txt
@@ -1,8 +1,13 @@
 pytest
 pytest-cov
 pytest-mock
+pytest-xdist
 coverage
 coveralls
 psutil
 netCDF4==1.5; platform_system == 'Windows' and python_version == "3.7"
 netCDF4==1.6; platform_system != 'Windows' or python_version != "3.7"
+dask[array]
+dask[diagnostics]
+dask[distributed]
+importlib-metadata==2.0; python_version == "3.7"
diff --git a/tests/test-models b/tests/test-models
index a19097ce..65b6e112 160000
--- a/tests/test-models
+++ b/tests/test-models
@@ -1,1 +1,1 @@
-Subproject commit a19097cee55a631652feeb5ecf74adcadea549e6
+Subproject commit 65b6e112e77ecb526ee1158791738c1dbe4d4085