Combining different fields and indices fails
If one combines 3D and 2D fields and makes use of indices, only the surface velocity currents are used to advect particles.  A quick fix is to add a line indices = indices.copy() in field.py:

diff --git a/parcels/field.py b/parcels/field.py
index 163c325..bb67c86 100644
--- a/parcels/field.py
+++ b/parcels/field.py
@@ -273,19 +273,32 @@ class Field(object):
                full_load is however sometimes necessary for plotting the fields.
         """
 
+        indices = indices.copy()
         if not isinstance(filenames, Iterable) or isinstance(filenames, str):
             filenames = [filenames]
         dimension_filename = dimension_filename if dimension_filename else filenames[0]
         if indices is None:
             indices = {}
         with NetcdfFileBuffer(dimension_filename, dimensions, indices) as filebuffer:

diff --git a/parcels/field.py b/parcels/field.py
index f49b351a..fa1851e0 100644
--- a/parcels/field.py
+++ b/parcels/field.py
@@ -152,8 +152,7 @@ class Field(object):
         if not isinstance(filenames, Iterable) or isinstance(filenames, str):
             filenames = [filenames]
         dimension_filename = dimension_filename if dimension_filename else filenames[0]
-        if indices is None:
-            indices = {}
+        indices = {} if indices is None else indices.copy()
         with NetcdfFileBuffer(dimension_filename, dimensions, indices) as filebuffer:
             lon, lat = filebuffer.read_lonlat
             depth = filebuffer.read_depth
diff --git a/tests/test_fieldset.py b/tests/test_fieldset.py
index c7ffa424..253cf9f9 100644
--- a/tests/test_fieldset.py
+++ b/tests/test_fieldset.py
@@ -92,7 +92,9 @@ def test_fieldset_from_file_subsets(indslon, indslat, tmpdir, filename='test_sub
     fieldsetfull = FieldSet.from_data(data, dimensions)
     fieldsetfull.write(filepath)
     indices = {'lon': indslon, 'lat': indslat}
+    indices_back = indices.copy()
     fieldsetsub = FieldSet.from_parcels(filepath, indices=indices)
+    assert indices == indices_back
     assert np.allclose(fieldsetsub.U.lon, fieldsetfull.U.grid.lon[indices['lon']])
     assert np.allclose(fieldsetsub.U.lat, fieldsetfull.U.grid.lat[indices['lat']])
     assert np.allclose(fieldsetsub.V.lon, fieldsetfull.V.grid.lon[indices['lon']])
['parcels/field.py', 'tests/test_fieldset.py']