diff --git a/parcels/codegenerator.py b/parcels/codegenerator.py
index 14cd6f797e..3b38e0f172 100644
--- a/parcels/codegenerator.py
+++ b/parcels/codegenerator.py
@@ -396,7 +396,7 @@ def generate(self, py_ast, funcvars):
                 funcvars.remove(kvar)
         self.ccode.body.insert(0, c.Value('ErrorCode', 'err'))
         if len(funcvars) > 0:
-            self.ccode.body.insert(0, c.Value("float", ", ".join(funcvars)))
+            self.ccode.body.insert(0, c.Value("type_coord", ", ".join(funcvars)))
         if len(transformer.tmp_vars) > 0:
             self.ccode.body.insert(0, c.Value("float", ", ".join(transformer.tmp_vars)))
 
diff --git a/parcels/include/index_search.h b/parcels/include/index_search.h
index d5a48181b2..d08c0ec3f8 100644
--- a/parcels/include/index_search.h
+++ b/parcels/include/index_search.h
@@ -10,6 +10,12 @@ extern "C" {
 
 #define CHECKERROR(res) do {if (res != SUCCESS) return res;} while (0)
 
+#ifdef DOUBLE_COORD_VARIABLES
+typedef double type_coord;
+#else
+typedef float type_coord;
+#endif
+
 typedef struct
 {
   int gtype;
@@ -38,7 +44,7 @@ typedef enum
     RECTILINEAR_Z_GRID=0, RECTILINEAR_S_GRID=1, CURVILINEAR_Z_GRID=2, CURVILINEAR_S_GRID=3
   } GridCode;
 
-static inline ErrorCode search_indices_vertical_z(float z, int zdim, float *zvals, int *zi, double *zeta)
+static inline ErrorCode search_indices_vertical_z(type_coord z, int zdim, float *zvals, int *zi, double *zeta)
 {
   if (z < zvals[0] || z > zvals[zdim-1]) {return ERROR_OUT_OF_BOUNDS;}
   while (*zi < zdim-1 && z > zvals[*zi+1]) ++(*zi);
@@ -49,7 +55,7 @@ static inline ErrorCode search_indices_vertical_z(float z, int zdim, float *zval
   return SUCCESS;
 }
 
-static inline ErrorCode search_indices_vertical_s(float z, int xdim, int ydim, int zdim, float *zvals,
+static inline ErrorCode search_indices_vertical_s(type_coord z, int xdim, int ydim, int zdim, float *zvals,
                                     int xi, int yi, int *zi, double xsi, double eta, double *zeta,
                                     int z4d, int ti, int tdim, double time, double t0, double t1)
 {
@@ -120,7 +126,7 @@ static inline void reconnect_bnd_indices(int *xi, int *yi, int xdim, int ydim, i
 }
 
 
-static inline ErrorCode search_indices_rectilinear(float x, float y, float z, CStructuredGrid *grid, GridCode gcode,
+static inline ErrorCode search_indices_rectilinear(type_coord x, type_coord y, type_coord z, CStructuredGrid *grid, GridCode gcode,
                                                    int *xi, int *yi, int *zi, double *xsi, double *eta, double *zeta,
                                                    int ti, double time, double t0, double t1)
 {
@@ -213,7 +219,7 @@ static inline ErrorCode search_indices_rectilinear(float x, float y, float z, CS
 }
 
 
-static inline ErrorCode search_indices_curvilinear(float x, float y, float z, CStructuredGrid *grid, GridCode gcode,
+static inline ErrorCode search_indices_curvilinear(type_coord x, type_coord y, type_coord z, CStructuredGrid *grid, GridCode gcode,
                                                    int *xi, int *yi, int *zi, double *xsi, double *eta, double *zeta,
                                                    int ti, double time, double t0, double t1)
 {
@@ -351,7 +357,7 @@ static inline ErrorCode search_indices_curvilinear(float x, float y, float z, CS
 /* Local linear search to update grid index
  * params ti, sizeT, time. t0, t1 are only used for 4D S grids
  * */
-static inline ErrorCode search_indices(float x, float y, float z, CStructuredGrid *grid,
+static inline ErrorCode search_indices(type_coord x, type_coord y, type_coord z, CStructuredGrid *grid,
                                        int *xi, int *yi, int *zi, double *xsi, double *eta, double *zeta,
                                        GridCode gcode, int ti, double time, double t0, double t1)
 {
diff --git a/parcels/include/parcels.h b/parcels/include/parcels.h
index 43fd5def8b..a80925fd3b 100644
--- a/parcels/include/parcels.h
+++ b/parcels/include/parcels.h
@@ -101,7 +101,7 @@ static inline ErrorCode spatial_interpolation_tracer_c_grid_3D(int xi, int yi, i
 }
 
 /* Linear interpolation along the time axis */
-static inline ErrorCode temporal_interpolation_structured_grid(float x, float y, float z, double time, CField *f, 
+static inline ErrorCode temporal_interpolation_structured_grid(type_coord x, type_coord y, type_coord z, double time, CField *f,
                                                                GridCode gcode, int *xi, int *yi, int *zi, int *ti,
                                                                float *value, int interp_method)
 {
@@ -303,7 +303,7 @@ static inline ErrorCode spatial_interpolation_UV_c_grid(double xsi, double eta,
 
 
 
-static inline ErrorCode temporal_interpolationUV_c_grid(float x, float y, float z, double time, CField *U, CField *V,
+static inline ErrorCode temporal_interpolationUV_c_grid(type_coord x, type_coord y, type_coord z, double time, CField *U, CField *V,
                                                          GridCode gcode, int *xi, int *yi, int *zi, int *ti,
                                                          float *u, float *v)
 {
@@ -490,7 +490,7 @@ static inline ErrorCode spatial_interpolation_UVW_c_grid(double xsi, double eta,
   return SUCCESS;
 }
 
-static inline ErrorCode temporal_interpolationUVW_c_grid(float x, float y, float z, double time, CField *U, CField *V, CField *W,
+static inline ErrorCode temporal_interpolationUVW_c_grid(type_coord x, type_coord y, type_coord z, double time, CField *U, CField *V, CField *W,
                                                          GridCode gcode, int *xi, int *yi, int *zi, int *ti,
                                                          float *u, float *v, float *w)
 {
@@ -540,7 +540,7 @@ static inline ErrorCode temporal_interpolationUVW_c_grid(float x, float y, float
 }
 
 
-static inline ErrorCode temporal_interpolation(float x, float y, float z, double time, CField *f, 
+static inline ErrorCode temporal_interpolation(type_coord x, type_coord y, type_coord z, double time, CField *f,
                                                void *vxi, void *vyi, void *vzi, void *vti,
                                                float *value, int interp_method)
 {
@@ -559,7 +559,7 @@ static inline ErrorCode temporal_interpolation(float x, float y, float z, double
   }
 }
 
-static inline ErrorCode temporal_interpolationUV(float x, float y, float z, double time,
+static inline ErrorCode temporal_interpolationUV(type_coord x, type_coord y, type_coord z, double time,
                                                  CField *U, CField *V,
                                                  void *vxi, void *vyi, void *vzi, void *vti,
                                                  float *valueU, float *valueV, int interp_method)
@@ -582,7 +582,7 @@ static inline ErrorCode temporal_interpolationUV(float x, float y, float z, doub
   }
 }
 
-static inline ErrorCode temporal_interpolationUVW(float x, float y, float z, double time,
+static inline ErrorCode temporal_interpolationUVW(type_coord x, type_coord y, type_coord z, double time,
                                                   CField *U, CField *V, CField *W,
                                                   void *vxi, void *vyi, void *vzi, void *vti,
                                                   float *valueU, float *valueV, float *valueW, int interp_method)
diff --git a/parcels/kernels/advection.py b/parcels/kernels/advection.py
index 0a3278f0cb..a1b1f2cbeb 100644
--- a/parcels/kernels/advection.py
+++ b/parcels/kernels/advection.py
@@ -57,8 +57,8 @@ def AdvectionRK45(particle, fieldset, time):
 
     Times-step dt is halved if error is larger than tolerance, and doubled
     if error is smaller than 1/10th of tolerance, with tolerance set to
-    1e-9 * dt by default."""
-    tol = [1e-9]
+    1e-5 * dt by default."""
+    tol = [1e-5]
     c = [1./4., 3./8., 12./13., 1., 1./2.]
     A = [[1./4., 0., 0., 0., 0.],
          [3./32., 9./32., 0., 0., 0.],
diff --git a/parcels/particle.py b/parcels/particle.py
index 4fd61f0e94..4bef1445c1 100644
--- a/parcels/particle.py
+++ b/parcels/particle.py
@@ -194,6 +194,12 @@ def __repr__(self):
     def delete(self):
         self.state = ErrorCode.Delete
 
+    @classmethod
+    def set_lonlatdepth_dtype(cls, dtype):
+        cls.lon.dtype = dtype
+        cls.lat.dtype = dtype
+        cls.depth.dtype = dtype
+
 
 class JITParticle(ScipyParticle):
     """Particle class for JIT-based (Just-In-Time) Particle objects
diff --git a/parcels/particleset.py b/parcels/particleset.py
index a4c1f03264..4ae3eb2f0f 100644
--- a/parcels/particleset.py
+++ b/parcels/particleset.py
@@ -5,6 +5,7 @@
 from parcels.particlefile import ParticleFile
 from parcels.tools.loggers import logger
 from parcels.grid import GridCode
+from parcels.field import NestedField, SummedField
 import numpy as np
 import progressbar
 import time as time_module
@@ -28,10 +29,13 @@ class ParticleSet(object):
     :param depth: Optional list of initial depth values for particles. Default is 0m
     :param time: Optional list of initial time values for particles. Default is fieldset.U.grid.time[0]
     :param repeatdt: Optional interval (in seconds) on which to repeat the release of the ParticleSet
+    :param lonlatdepth_dtype: Floating precision for lon, lat, depth particle coordinates.
+           It is either np.float32 or np.float64. Default is np.float32 if fieldset.U.interp_method is 'linear'
+           and np.float64 if the interpolation method is 'cgrid_velocity'
     Other Variables can be initialised using further arguments (e.g. v=... for a Variable named 'v')
     """
 
-    def __init__(self, fieldset, pclass=JITParticle, lon=None, lat=None, depth=None, time=None, repeatdt=None, **kwargs):
+    def __init__(self, fieldset, pclass=JITParticle, lon=None, lat=None, depth=None, time=None, repeatdt=None, lonlatdepth_dtype=None, **kwargs):
         self.fieldset = fieldset
         self.fieldset.check_complete()
 
@@ -82,6 +86,14 @@ def convert_to_list(var):
             self.repeatpclass = pclass
             self.repeatkwargs = kwargs
 
+        if lonlatdepth_dtype is None:
+            self.lonlatdepth_dtype = self.lonlatdepth_dtype_from_field_interp_method(fieldset.U)
+        else:
+            self.lonlatdepth_dtype = lonlatdepth_dtype
+        assert self.lonlatdepth_dtype in [np.float32, np.float64], \
+            'lon lat depth precision should be set to either np.float32 or np.float64'
+        JITParticle.set_lonlatdepth_dtype(self.lonlatdepth_dtype)
+
         size = len(lon)
         self.particles = np.empty(size, dtype=pclass)
         self.ptype = pclass.getPType()
@@ -143,8 +155,10 @@ def from_line(cls, fieldset, pclass, start, finish, size, depth=None, time=None,
         :param time: Optional start time value for particles. Default is fieldset.U.time[0]
         :param repeatdt: Optional interval (in seconds) on which to repeat the release of the ParticleSet
         """
-        lon = np.linspace(start[0], finish[0], size, dtype=np.float32)
-        lat = np.linspace(start[1], finish[1], size, dtype=np.float32)
+
+        lonlat_type = cls.lonlatdepth_dtype_from_field_interp_method(fieldset.U)
+        lon = np.linspace(start[0], finish[0], size, dtype=lonlat_type)
+        lat = np.linspace(start[1], finish[1], size, dtype=lonlat_type)
         if type(depth) in [int, float]:
             depth = [depth] * size
         return cls(fieldset=fieldset, pclass=pclass, lon=lon, lat=lat, depth=depth, time=time, repeatdt=repeatdt)
@@ -201,6 +215,17 @@ def from_field(cls, fieldset, pclass, start_field, size, mode='monte_carlo', dep
 
         return cls(fieldset=fieldset, pclass=pclass, lon=lon, lat=lat, depth=depth, time=time, repeatdt=repeatdt)
 
+    @staticmethod
+    def lonlatdepth_dtype_from_field_interp_method(field):
+        if type(field) in [SummedField, NestedField]:
+            for f in field:
+                if f.interp_method == 'cgrid_velocity':
+                    return np.float64
+        else:
+            if field.interp_method == 'cgrid_velocity':
+                return np.float64
+        return np.float32
+
     @property
     def size(self):
         return self.particles.size
@@ -290,7 +315,8 @@ def execute(self, pyfunc=AdvectionRK4, endtime=None, runtime=None, dt=1.,
             # Prepare JIT kernel execution
             if self.ptype.uses_jit:
                 self.kernel.remove_lib()
-                self.kernel.compile(compiler=GNUCompiler())
+                cppargs = ['-DDOUBLE_COORD_VARIABLES'] if self.lonlatdepth_dtype == np.float64 else None
+                self.kernel.compile(compiler=GNUCompiler(cppargs=cppargs))
                 self.kernel.load_lib()
 
         # Convert all time variables to seconds
diff --git a/tests/customed_header.h b/tests/customed_header.h
index e861c36c49..a9b2072bf3 100644
--- a/tests/customed_header.h
+++ b/tests/customed_header.h
@@ -1,5 +1,5 @@
 
-static inline void func(CField *f, float *lon, float *dt)
+static inline void func(CField *f, double *lon, float *dt)
 {
   float (*data)[f->xdim] = (float (*)[f->xdim]) f->data;
   float u = data[2][1];
diff --git a/tests/test_advection.py b/tests/test_advection.py
index 4897c902d4..216350f36f 100644
--- a/tests/test_advection.py
+++ b/tests/test_advection.py
@@ -47,17 +47,17 @@ def test_advection_zonal(lon, lat, depth, mode, npart=10):
     assert fieldset2D.U.creation_log == 'from_data'
 
     pset2D = ParticleSet(fieldset2D, pclass=ptype[mode],
-                         lon=np.zeros(npart, dtype=np.float32) + 20.,
-                         lat=np.linspace(0, 80, npart, dtype=np.float32))
+                         lon=np.zeros(npart) + 20.,
+                         lat=np.linspace(0, 80, npart))
     pset2D.execute(AdvectionRK4, runtime=delta(hours=2), dt=delta(seconds=30))
     assert (np.diff(np.array([p.lon for p in pset2D])) > 1.e-4).all()
 
     dimensions['depth'] = depth
     fieldset3D = FieldSet.from_data(data3D, dimensions, mesh='spherical', transpose=True)
     pset3D = ParticleSet(fieldset3D, pclass=ptype[mode],
-                         lon=np.zeros(npart, dtype=np.float32) + 20.,
-                         lat=np.linspace(0, 80, npart, dtype=np.float32),
-                         depth=np.zeros(npart, dtype=np.float32) + 10.)
+                         lon=np.zeros(npart) + 20.,
+                         lat=np.linspace(0, 80, npart),
+                         depth=np.zeros(npart) + 10.)
     pset3D.execute(AdvectionRK4, runtime=delta(hours=2), dt=delta(seconds=30))
     assert (np.diff(np.array([p.lon for p in pset3D])) > 1.e-4).all()
 
@@ -73,8 +73,8 @@ def test_advection_meridional(lon, lat, mode, npart=10):
     fieldset = FieldSet.from_data(data, dimensions, mesh='spherical', transpose=True)
 
     pset = ParticleSet(fieldset, pclass=ptype[mode],
-                       lon=np.linspace(-60, 60, npart, dtype=np.float32),
-                       lat=np.linspace(0, 30, npart, dtype=np.float32))
+                       lon=np.linspace(-60, 60, npart),
+                       lat=np.linspace(0, 30, npart))
     delta_lat = np.diff(np.array([p.lat for p in pset]))
     pset.execute(AdvectionRK4, runtime=delta(hours=2), dt=delta(seconds=30))
     assert np.allclose(np.diff(np.array([p.lat for p in pset])), delta_lat, rtol=1.e-4)
@@ -94,9 +94,9 @@ def test_advection_3D(mode, npart=11):
     fieldset = FieldSet.from_data(data, dimensions, mesh='flat', transpose=True)
 
     pset = ParticleSet(fieldset, pclass=ptype[mode],
-                       lon=np.zeros(npart, dtype=np.float32),
-                       lat=np.zeros(npart, dtype=np.float32) + 1e2,
-                       depth=np.linspace(0, 1, npart, dtype=np.float32))
+                       lon=np.zeros(npart),
+                       lat=np.zeros(npart) + 1e2,
+                       depth=np.linspace(0, 1, npart))
     time = delta(hours=2).total_seconds()
     pset.execute(AdvectionRK4, runtime=time, dt=delta(seconds=30))
     assert np.allclose([p.depth*time for p in pset], [p.lon for p in pset], atol=1.e-1)
@@ -182,8 +182,8 @@ def fieldset_stationary(xdim=100, ydim=100, maxtime=delta(hours=6)):
     ('RK45', 1e-5)])
 def test_stationary_eddy(fieldset_stationary, mode, method, rtol, npart=1):
     fieldset = fieldset_stationary
-    lon = np.linspace(12000, 21000, npart, dtype=np.float32)
-    lat = np.linspace(12500, 12500, npart, dtype=np.float32)
+    lon = np.linspace(12000, 21000, npart)
+    lat = np.linspace(12500, 12500, npart)
     pset = ParticleSet(fieldset, pclass=ptype[mode], lon=lon, lat=lat)
     endtime = delta(hours=6).total_seconds()
     pset.execute(kernel[method], dt=delta(minutes=3), endtime=endtime)
@@ -195,9 +195,9 @@ def test_stationary_eddy(fieldset_stationary, mode, method, rtol, npart=1):
 
 @pytest.mark.parametrize('mode', ['scipy', 'jit'])
 def test_stationary_eddy_vertical(mode, npart=1):
-    lon = np.linspace(12000, 21000, npart, dtype=np.float32)
-    lat = np.linspace(10000, 20000, npart, dtype=np.float32)
-    depth = np.linspace(12500, 12500, npart, dtype=np.float32)
+    lon = np.linspace(12000, 21000, npart)
+    lat = np.linspace(10000, 20000, npart)
+    depth = np.linspace(12500, 12500, npart)
     endtime = delta(hours=6).total_seconds()
 
     xdim = ydim = 100
@@ -263,8 +263,8 @@ def fieldset_moving(xdim=100, ydim=100, maxtime=delta(hours=6)):
     ('RK45', 1e-5)])
 def test_moving_eddy(fieldset_moving, mode, method, rtol, npart=1):
     fieldset = fieldset_moving
-    lon = np.linspace(12000, 21000, npart, dtype=np.float32)
-    lat = np.linspace(12500, 12500, npart, dtype=np.float32)
+    lon = np.linspace(12000, 21000, npart)
+    lat = np.linspace(12500, 12500, npart)
     pset = ParticleSet(fieldset, pclass=ptype[mode], lon=lon, lat=lat)
     endtime = delta(hours=6).total_seconds()
     pset.execute(kernel[method], dt=delta(minutes=3), endtime=endtime)
@@ -307,8 +307,8 @@ def fieldset_decaying(xdim=100, ydim=100, maxtime=delta(hours=6)):
     ('RK45', 1e-5)])
 def test_decaying_eddy(fieldset_decaying, mode, method, rtol, npart=1):
     fieldset = fieldset_decaying
-    lon = np.linspace(12000, 21000, npart, dtype=np.float32)
-    lat = np.linspace(12500, 12500, npart, dtype=np.float32)
+    lon = np.linspace(12000, 21000, npart)
+    lat = np.linspace(12500, 12500, npart)
     pset = ParticleSet(fieldset, pclass=ptype[mode], lon=lon, lat=lat)
     endtime = delta(hours=6).total_seconds()
     pset.execute(kernel[method], dt=delta(minutes=3), endtime=endtime)
@@ -344,8 +344,8 @@ def test_decaying_eddy(fieldset_decaying, mode, method, rtol, npart=1):
 
     npart = args.particles
     pset = ParticleSet(fieldset, pclass=ptype[args.mode],
-                       lon=np.linspace(4000, 21000, npart, dtype=np.float32),
-                       lat=np.linspace(12500, 12500, npart, dtype=np.float32))
+                       lon=np.linspace(4000, 21000, npart),
+                       lat=np.linspace(12500, 12500, npart))
     if args.verbose:
         print("Initial particle positions:\n%s" % pset)
     pset.execute(kernel[args.method], dt=delta(minutes=3), runtime=delta(hours=6))
diff --git a/tests/test_fieldset.py b/tests/test_fieldset.py
index d84f5f5283..a595d9ad8d 100644
--- a/tests/test_fieldset.py
+++ b/tests/test_fieldset.py
@@ -436,7 +436,7 @@ def generate_dataset(xdim, ydim, zdim=1, tdim=1):
     pset = ParticleSet(fieldset, JITParticle, 0, 0)
 
     pset.execute(AdvectionRK4, dt=1)
-    assert pset[0].lon == 4.5 and pset[0].lat == 10
+    assert np.allclose(pset[0].lon, 4.5) and np.allclose(pset[0].lat, 10)
 
 
 def test_fieldset_from_data_gridtypes(xdim=20, ydim=10, zdim=4):
diff --git a/tests/test_fieldset_sampling.py b/tests/test_fieldset_sampling.py
index 415e2d026e..190c6ad54a 100644
--- a/tests/test_fieldset_sampling.py
+++ b/tests/test_fieldset_sampling.py
@@ -213,13 +213,13 @@ def test_fieldset_sample_particle(mode, k_sample_uv, lat_flip, npart=120):
 
     fieldset = FieldSet.from_data(data, dimensions, mesh='flat', transpose=True)
 
-    lon = np.linspace(-170, 170, npart, dtype=np.float32)
-    lat = np.linspace(-80, 80, npart, dtype=np.float32)
-    pset = ParticleSet(fieldset, pclass=pclass(mode), lon=lon, lat=np.zeros(npart, dtype=np.float32) + 70.)
+    lon = np.linspace(-170, 170, npart)
+    lat = np.linspace(-80, 80, npart)
+    pset = ParticleSet(fieldset, pclass=pclass(mode), lon=lon, lat=np.zeros(npart) + 70.)
     pset.execute(pset.Kernel(k_sample_uv), endtime=1., dt=1.)
     assert np.allclose(np.array([p.v for p in pset]), lon, rtol=1e-6)
 
-    pset = ParticleSet(fieldset, pclass=pclass(mode), lat=lat, lon=np.zeros(npart, dtype=np.float32) - 45.)
+    pset = ParticleSet(fieldset, pclass=pclass(mode), lat=lat, lon=np.zeros(npart) - 45.)
     pset.execute(pset.Kernel(k_sample_uv), endtime=1., dt=1.)
     assert np.allclose(np.array([p.u for p in pset]), lat, rtol=1e-6)
 
@@ -228,14 +228,14 @@ def test_fieldset_sample_particle(mode, k_sample_uv, lat_flip, npart=120):
 def test_fieldset_sample_geographic(fieldset_geometric, mode, k_sample_uv, npart=120):
     """ Sample a fieldset with conversion to geographic units (degrees). """
     fieldset = fieldset_geometric
-    lon = np.linspace(-170, 170, npart, dtype=np.float32)
-    lat = np.linspace(-80, 80, npart, dtype=np.float32)
+    lon = np.linspace(-170, 170, npart)
+    lat = np.linspace(-80, 80, npart)
 
-    pset = ParticleSet(fieldset, pclass=pclass(mode), lon=lon, lat=np.zeros(npart, dtype=np.float32) + 70.)
+    pset = ParticleSet(fieldset, pclass=pclass(mode), lon=lon, lat=np.zeros(npart) + 70.)
     pset.execute(pset.Kernel(k_sample_uv), endtime=1., dt=1.)
     assert np.allclose(np.array([p.v for p in pset]), lon, rtol=1e-6)
 
-    pset = ParticleSet(fieldset, pclass=pclass(mode), lat=lat, lon=np.zeros(npart, dtype=np.float32) - 45.)
+    pset = ParticleSet(fieldset, pclass=pclass(mode), lat=lat, lon=np.zeros(npart) - 45.)
     pset.execute(pset.Kernel(k_sample_uv), endtime=1., dt=1.)
     assert np.allclose(np.array([p.u for p in pset]), lat, rtol=1e-6)
 
@@ -244,14 +244,14 @@ def test_fieldset_sample_geographic(fieldset_geometric, mode, k_sample_uv, npart
 def test_fieldset_sample_geographic_polar(fieldset_geometric_polar, mode, k_sample_uv, npart=120):
     """ Sample a fieldset with conversion to geographic units and a pole correction. """
     fieldset = fieldset_geometric_polar
-    lon = np.linspace(-170, 170, npart, dtype=np.float32)
-    lat = np.linspace(-80, 80, npart, dtype=np.float32)
+    lon = np.linspace(-170, 170, npart)
+    lat = np.linspace(-80, 80, npart)
 
-    pset = ParticleSet(fieldset, pclass=pclass(mode), lon=lon, lat=np.zeros(npart, dtype=np.float32) + 70.)
+    pset = ParticleSet(fieldset, pclass=pclass(mode), lon=lon, lat=np.zeros(npart) + 70.)
     pset.execute(pset.Kernel(k_sample_uv), endtime=1., dt=1.)
     assert np.allclose(np.array([p.v for p in pset]), lon, rtol=1e-6)
 
-    pset = ParticleSet(fieldset, pclass=pclass(mode), lat=lat, lon=np.zeros(npart, dtype=np.float32) - 45.)
+    pset = ParticleSet(fieldset, pclass=pclass(mode), lat=lat, lon=np.zeros(npart) - 45.)
     pset.execute(pset.Kernel(k_sample_uv), endtime=1., dt=1.)
     # Note: 1.e-2 is a very low rtol, so there seems to be a rather
     # large sampling error for the JIT correction.
diff --git a/tests/test_kernel_execution.py b/tests/test_kernel_execution.py
index f586fd217e..5081e3ee3b 100644
--- a/tests/test_kernel_execution.py
+++ b/tests/test_kernel_execution.py
@@ -35,8 +35,8 @@ def fieldset(xdim=20, ydim=20):
 ])
 def test_execution_endtime(fieldset, mode, start, end, substeps, dt, npart=10):
     pset = ParticleSet(fieldset, pclass=ptype[mode], time=start,
-                       lon=np.linspace(0, 1, npart, dtype=np.float32),
-                       lat=np.linspace(1, 0, npart, dtype=np.float32))
+                       lon=np.linspace(0, 1, npart),
+                       lat=np.linspace(1, 0, npart))
     pset.execute(DoNothing, endtime=end, dt=dt)
     assert np.allclose(np.array([p.time for p in pset]), end)
 
@@ -52,8 +52,8 @@ def test_execution_endtime(fieldset, mode, start, end, substeps, dt, npart=10):
 ])
 def test_execution_runtime(fieldset, mode, start, end, substeps, dt, npart=10):
     pset = ParticleSet(fieldset, pclass=ptype[mode], time=start,
-                       lon=np.linspace(0, 1, npart, dtype=np.float32),
-                       lat=np.linspace(1, 0, npart, dtype=np.float32))
+                       lon=np.linspace(0, 1, npart),
+                       lat=np.linspace(1, 0, npart))
     t_step = abs(end - start) / substeps
     for _ in range(substeps):
         pset.execute(DoNothing, runtime=t_step, dt=dt)
@@ -66,8 +66,8 @@ def test_execution_runtime(fieldset, mode, start, end, substeps, dt, npart=10):
 def test_pset_execute_dt_0(fieldset, mode, time, dt, npart=2):
     def SetLat(particle, fieldset, time):
         particle.lat = .6
-    lon = np.linspace(0, 1, npart, dtype=np.float32)
-    lat = np.linspace(1, 0, npart, dtype=np.float32)
+    lon = np.linspace(0, 1, npart)
+    lat = np.linspace(1, 0, npart)
 
     pset = ParticleSet(fieldset, pclass=ptype[mode], lon=lon, lat=lat)
     pset.execute(SetLat, endtime=time, dt=dt)
@@ -91,8 +91,8 @@ def TimedFail(particle, fieldset, time):
             return ErrorCode.Success
 
     pset = ParticleSet(fieldset, pclass=ptype[mode],
-                       lon=np.linspace(0, 1, npart, dtype=np.float32),
-                       lat=np.linspace(1, 0, npart, dtype=np.float32))
+                       lon=np.linspace(0, 1, npart),
+                       lat=np.linspace(1, 0, npart))
     error_thrown = False
     try:
         pset.execute(TimedFail, endtime=20., dt=2.)
@@ -112,8 +112,8 @@ def PythonFail(particle, fieldset, time):
             return ErrorCode.Success
 
     pset = ParticleSet(fieldset, pclass=ptype[mode],
-                       lon=np.linspace(0, 1, npart, dtype=np.float32),
-                       lat=np.linspace(1, 0, npart, dtype=np.float32))
+                       lon=np.linspace(0, 1, npart),
+                       lat=np.linspace(1, 0, npart))
     error_thrown = False
     try:
         pset.execute(PythonFail, endtime=20., dt=2.)
@@ -131,8 +131,8 @@ def MoveRight(particle, fieldset, time):
         particle.lon += 0.1
 
     pset = ParticleSet(fieldset, pclass=ptype[mode],
-                       lon=np.linspace(0, 1, npart, dtype=np.float32),
-                       lat=np.linspace(1, 0, npart, dtype=np.float32))
+                       lon=np.linspace(0, 1, npart),
+                       lat=np.linspace(1, 0, npart))
     error_thrown = False
     try:
         pset.execute(MoveRight, endtime=10., dt=1.)
@@ -152,8 +152,8 @@ def MoveRight(particle, fieldset, time):
     def MoveLeft(particle, fieldset, time):
         particle.lon -= 1.
 
-    lon = np.linspace(0.05, 0.95, npart, dtype=np.float32)
-    lat = np.linspace(1, 0, npart, dtype=np.float32)
+    lon = np.linspace(0.05, 0.95, npart)
+    lat = np.linspace(1, 0, npart)
     pset = ParticleSet(fieldset, pclass=ptype[mode], lon=lon, lat=lat)
     pset.execute(MoveRight, endtime=10., dt=1.,
                  recovery={ErrorCode.ErrorOutOfBounds: MoveLeft})
@@ -171,8 +171,8 @@ def MoveRight(particle, fieldset, time):
     def DeleteMe(particle, fieldset, time):
         particle.delete()
 
-    lon = np.linspace(0.05, 0.95, npart, dtype=np.float32)
-    lat = np.linspace(1, 0, npart, dtype=np.float32)
+    lon = np.linspace(0.05, 0.95, npart)
+    lat = np.linspace(1, 0, npart)
     pset = ParticleSet(fieldset, pclass=ptype[mode], lon=lon, lat=lat)
     pset.execute(MoveRight, endtime=10., dt=1.,
                  recovery={ErrorCode.ErrorOutOfBounds: DeleteMe})
diff --git a/tests/test_kernel_language.py b/tests/test_kernel_language.py
index 874bc9816e..b09e851c4d 100644
--- a/tests/test_kernel_language.py
+++ b/tests/test_kernel_language.py
@@ -41,8 +41,8 @@ def test_expression_int(fieldset, mode, name, expr, result, npart=10):
     class TestParticle(ptype[mode]):
         p = Variable('p', dtype=np.float32)
     pset = ParticleSet(fieldset, pclass=TestParticle,
-                       lon=np.linspace(0., 1., npart, dtype=np.float32),
-                       lat=np.zeros(npart, dtype=np.float32) + 0.5)
+                       lon=np.linspace(0., 1., npart),
+                       lat=np.zeros(npart) + 0.5)
     pset.execute(expr_kernel('Test%s' % name, pset, expr), endtime=1., dt=1.)
     assert(np.array([result == particle.p for particle in pset]).all())
 
@@ -60,8 +60,8 @@ def test_expression_float(fieldset, mode, name, expr, result, npart=10):
     class TestParticle(ptype[mode]):
         p = Variable('p', dtype=np.float32)
     pset = ParticleSet(fieldset, pclass=TestParticle,
-                       lon=np.linspace(0., 1., npart, dtype=np.float32),
-                       lat=np.zeros(npart, dtype=np.float32) + 0.5)
+                       lon=np.linspace(0., 1., npart),
+                       lat=np.zeros(npart) + 0.5)
     pset.execute(expr_kernel('Test%s' % name, pset, expr), endtime=1., dt=1.)
     assert(np.array([result == particle.p for particle in pset]).all())
 
@@ -83,8 +83,8 @@ def test_expression_bool(fieldset, mode, name, expr, result, npart=10):
     class TestParticle(ptype[mode]):
         p = Variable('p', dtype=np.float32)
     pset = ParticleSet(fieldset, pclass=TestParticle,
-                       lon=np.linspace(0., 1., npart, dtype=np.float32),
-                       lat=np.zeros(npart, dtype=np.float32) + 0.5)
+                       lon=np.linspace(0., 1., npart),
+                       lat=np.zeros(npart) + 0.5)
     pset.execute(expr_kernel('Test%s' % name, pset, expr), endtime=1., dt=1.)
     if mode == 'jit':
         assert(np.array([result == (particle.p == 1) for particle in pset]).all())
@@ -221,8 +221,8 @@ def test_random_float(fieldset, mode, rngfunc, rngargs, npart=10):
     class TestParticle(ptype[mode]):
         p = Variable('p', dtype=np.float32 if rngfunc == 'randint' else np.float32)
     pset = ParticleSet(fieldset, pclass=TestParticle,
-                       lon=np.linspace(0., 1., npart, dtype=np.float32),
-                       lat=np.zeros(npart, dtype=np.float32) + 0.5)
+                       lon=np.linspace(0., 1., npart),
+                       lat=np.zeros(npart) + 0.5)
     series = random_series(npart, rngfunc, rngargs, mode)
     kernel = expr_kernel('TestRandom_%s' % rngfunc, pset,
                          'random.%s(%s)' % (rngfunc, ', '.join([str(a) for a in rngargs])))
@@ -233,7 +233,9 @@ class TestParticle(ptype[mode]):
 @pytest.mark.parametrize('mode', ['scipy', 'jit'])
 @pytest.mark.parametrize('c_inc', ['str', 'file'])
 def test_c_kernel(fieldset, mode, c_inc):
-    pset = ParticleSet(fieldset, pclass=ptype[mode], lon=[0.5], lat=[0])
+    coord_type = np.float32 if c_inc == 'str' else np.float64
+    pset = ParticleSet(fieldset, pclass=ptype[mode], lon=[0.5], lat=[0],
+                       lonlatdepth_dtype=coord_type)
 
     def func(U, lon, dt):
         u = U.data[0, 2, 1]
diff --git a/tests/test_particle_sets.py b/tests/test_particle_sets.py
index 3ad1fa2ed5..ef2d7b88dc 100644
--- a/tests/test_particle_sets.py
+++ b/tests/test_particle_sets.py
@@ -100,8 +100,8 @@ def test_pset_create_field_curvi(npart=100):
 
 @pytest.mark.parametrize('mode', ['scipy', 'jit'])
 def test_pset_create_with_time(fieldset, mode, npart=100):
-    lon = np.linspace(0, 1, npart, dtype=np.float32)
-    lat = np.linspace(1, 0, npart, dtype=np.float32)
+    lon = np.linspace(0, 1, npart)
+    lat = np.linspace(1, 0, npart)
     time = 5.
     pset = ParticleSet(fieldset, lon=lon, lat=lat, pclass=ptype[mode], time=time)
     assert np.allclose([p.time for p in pset], time, rtol=1e-12)
@@ -208,8 +208,8 @@ def __init__(self, *args, **kwargs):
             self.n = 2
 
     pset = ParticleSet(fieldset, pclass=TestParticle,
-                       lon=np.linspace(0, 1, npart, dtype=np.float32),
-                       lat=np.linspace(1, 0, npart, dtype=np.float32))
+                       lon=np.linspace(0, 1, npart),
+                       lat=np.linspace(1, 0, npart))
     assert(pset.size == 100)
     # FIXME: The float test fails with a conversion error of 1.e-8
     # assert np.allclose([p.p - 0.33 for p in pset], np.zeros(npart), rtol=1e-12)
@@ -218,9 +218,9 @@ def __init__(self, *args, **kwargs):
 
 @pytest.mark.parametrize('mode', ['scipy', 'jit'])
 def test_pset_add_explicit(fieldset, mode, npart=100):
-    lon = np.linspace(0, 1, npart, dtype=np.float32)
-    lat = np.linspace(1, 0, npart, dtype=np.float32)
-    pset = ParticleSet(fieldset, lon=[], lat=[], pclass=ptype[mode])
+    lon = np.linspace(0, 1, npart)
+    lat = np.linspace(1, 0, npart)
+    pset = ParticleSet(fieldset, lon=[], lat=[], pclass=ptype[mode], lonlatdepth_dtype=np.float64)
     for i in range(npart):
         particle = ptype[mode](lon=lon[i], lat=lat[i], fieldset=fieldset)
         pset.add(particle)
@@ -257,11 +257,11 @@ def AddLat(particle, fieldset, time):
 @pytest.mark.parametrize('mode', ['scipy', 'jit'])
 def test_pset_merge_inplace(fieldset, mode, npart=100):
     pset1 = ParticleSet(fieldset, pclass=ptype[mode],
-                        lon=np.linspace(0, 1, npart, dtype=np.float32),
-                        lat=np.linspace(1, 0, npart, dtype=np.float32))
+                        lon=np.linspace(0, 1, npart),
+                        lat=np.linspace(1, 0, npart))
     pset2 = ParticleSet(fieldset, pclass=ptype[mode],
-                        lon=np.linspace(0, 1, npart, dtype=np.float32),
-                        lat=np.linspace(0, 1, npart, dtype=np.float32))
+                        lon=np.linspace(0, 1, npart),
+                        lat=np.linspace(0, 1, npart))
     assert(pset1.size == 100)
     assert(pset2.size == 100)
     pset1.add(pset2)
@@ -272,11 +272,11 @@ def test_pset_merge_inplace(fieldset, mode, npart=100):
 @pytest.mark.parametrize('mode', ['scipy', 'jit'])
 def test_pset_merge_duplicate(fieldset, mode, npart=100):
     pset1 = ParticleSet(fieldset, pclass=ptype[mode],
-                        lon=np.linspace(0, 1, npart, dtype=np.float32),
-                        lat=np.linspace(1, 0, npart, dtype=np.float32))
+                        lon=np.linspace(0, 1, npart),
+                        lat=np.linspace(1, 0, npart))
     pset2 = ParticleSet(fieldset, pclass=ptype[mode],
-                        lon=np.linspace(0, 1, npart, dtype=np.float32),
-                        lat=np.linspace(0, 1, npart, dtype=np.float32))
+                        lon=np.linspace(0, 1, npart),
+                        lat=np.linspace(0, 1, npart))
     pset3 = pset1 + pset2
     assert(pset1.size == 100)
     assert(pset2.size == 100)
@@ -285,9 +285,9 @@ def test_pset_merge_duplicate(fieldset, mode, npart=100):
 
 @pytest.mark.parametrize('mode', ['scipy', 'jit'])
 def test_pset_remove_index(fieldset, mode, npart=100):
-    lon = np.linspace(0, 1, npart, dtype=np.float32)
-    lat = np.linspace(1, 0, npart, dtype=np.float32)
-    pset = ParticleSet(fieldset, lon=lon, lat=lat, pclass=ptype[mode])
+    lon = np.linspace(0, 1, npart)
+    lat = np.linspace(1, 0, npart)
+    pset = ParticleSet(fieldset, lon=lon, lat=lat, pclass=ptype[mode], lonlatdepth_dtype=np.float64)
     for ilon, ilat in zip(lon[::-1], lat[::-1]):
         p = pset.remove(-1)
         assert(p.lon == ilon)
@@ -298,8 +298,8 @@ def test_pset_remove_index(fieldset, mode, npart=100):
 @pytest.mark.xfail(reason="Particle removal has not been implemented yet")
 @pytest.mark.parametrize('mode', ['scipy', 'jit'])
 def test_pset_remove_particle(fieldset, mode, npart=100):
-    lon = np.linspace(0, 1, npart, dtype=np.float32)
-    lat = np.linspace(1, 0, npart, dtype=np.float32)
+    lon = np.linspace(0, 1, npart)
+    lat = np.linspace(1, 0, npart)
     pset = ParticleSet(fieldset, lon=lon, lat=lat, pclass=ptype[mode])
     for ilon, ilat in zip(lon[::-1], lat[::-1]):
         p = pset.remove(pset[-1])
@@ -315,8 +315,8 @@ def DeleteKernel(particle, fieldset, time):
             particle.delete()
 
     pset = ParticleSet(fieldset, pclass=ptype[mode],
-                       lon=np.linspace(0, 1, npart, dtype=np.float32),
-                       lat=np.linspace(1, 0, npart, dtype=np.float32))
+                       lon=np.linspace(0, 1, npart),
+                       lat=np.linspace(1, 0, npart))
     pset.execute(pset.Kernel(DeleteKernel), endtime=1., dt=1.0)
     assert(pset.size == 40)
 
@@ -327,8 +327,8 @@ def AddLat(particle, fieldset, time):
         particle.lat += 0.1
 
     pset = ParticleSet(fieldset, pclass=ptype[mode],
-                       lon=np.linspace(0, 1, npart, dtype=np.float32),
-                       lat=np.zeros(npart, dtype=np.float32))
+                       lon=np.linspace(0, 1, npart),
+                       lat=np.zeros(npart))
     k_add = pset.Kernel(AddLat)
     for _ in range(n):
         pset.execute(k_add, runtime=1., dt=1.0)
@@ -341,8 +341,8 @@ def AddLat(particle, fieldset, time):
         particle.lat += 0.1
 
     pset = ParticleSet(fieldset, pclass=ptype[mode],
-                       lon=np.linspace(0, 1, npart, dtype=np.float32),
-                       lat=np.zeros(npart, dtype=np.float32))
+                       lon=np.linspace(0, 1, npart),
+                       lat=np.zeros(npart))
     k_add = pset.Kernel(AddLat)
     for _ in range(n):
         pset.execute(k_add, runtime=1., dt=1.0)
@@ -370,8 +370,8 @@ def test_density(fieldset, mode, area_scale):
 def test_pfile_array_remove_particles(fieldset, mode, tmpdir, npart=10):
     filepath = tmpdir.join("pfile_array_remove_particles")
     pset = ParticleSet(fieldset, pclass=ptype[mode],
-                       lon=np.linspace(0, 1, npart, dtype=np.float32),
-                       lat=0.5*np.ones(npart, dtype=np.float32))
+                       lon=np.linspace(0, 1, npart),
+                       lat=0.5*np.ones(npart))
     pfile = pset.ParticleFile(filepath)
     pfile.write(pset, 0)
     pset.remove(3)
@@ -383,8 +383,8 @@ def test_pfile_array_remove_all_particles(fieldset, mode, tmpdir, npart=10):
 
     filepath = tmpdir.join("pfile_array_remove_particles")
     pset = ParticleSet(fieldset, pclass=ptype[mode],
-                       lon=np.linspace(0, 1, npart, dtype=np.float32),
-                       lat=0.5*np.ones(npart, dtype=np.float32))
+                       lon=np.linspace(0, 1, npart),
+                       lat=0.5*np.ones(npart))
     pfile = pset.ParticleFile(filepath)
     pfile.write(pset, 0)
     for _ in range(npart):
@@ -405,8 +405,8 @@ def Update_v(particle, fieldset, time):
     class MyParticle(ptype[mode]):
         v_once = Variable('v_once', dtype=np.float32, initial=0., to_write='once')
         age = Variable('age', dtype=np.float32, initial=0.)
-    lon = np.linspace(0, 1, npart, dtype=np.float32)
-    lat = np.linspace(1, 0, npart, dtype=np.float32)
+    lon = np.linspace(0, 1, npart)
+    lat = np.linspace(1, 0, npart)
     pset = ParticleSet(fieldset, pclass=MyParticle, lon=lon, lat=lat, repeatdt=0.1)
     pset.execute(pset.Kernel(Update_v), endtime=1, dt=0.1,
                  output_file=pset.ParticleFile(name=filepath, outputdt=0.1))
@@ -428,8 +428,8 @@ def move_west(particle, fieldset, time):
     def DeleteP(particle, fieldset, time):
         particle.delete()
 
-    lon = np.linspace(0.05, 0.95, npart, dtype=np.float32)
-    lat = np.linspace(0.95, 0.05, npart, dtype=np.float32)
+    lon = np.linspace(0.05, 0.95, npart)
+    lat = np.linspace(0.95, 0.05, npart)
 
     (dt, runtime) = (0.1, 0.8)
     lon_end = lon - runtime/dt*0.1
diff --git a/tests/test_particles.py b/tests/test_particles.py
index 2bdb64fc5e..7155309fe1 100644
--- a/tests/test_particles.py
+++ b/tests/test_particles.py
@@ -24,8 +24,8 @@ class TestParticle(ptype[mode]):
         p_double = Variable('p_double', dtype=np.float64, initial=11.)
         p_int = Variable('p_int', dtype=np.int32, initial=12.)
     pset = ParticleSet(fieldset, pclass=TestParticle,
-                       lon=np.linspace(0, 1, npart, dtype=np.float32),
-                       lat=np.linspace(1, 0, npart, dtype=np.float32))
+                       lon=np.linspace(0, 1, npart),
+                       lat=np.linspace(1, 0, npart))
 
     def addOne(particle, fieldset, time):
         particle.p_float += 1.
@@ -65,25 +65,28 @@ class TestParticle(ptype[mode]):
 
 
 @pytest.mark.parametrize('mode', ['scipy', 'jit'])
-def test_variable_init_relative(fieldset, mode, npart=10):
+@pytest.mark.parametrize('coord_type', [np.float32, np.float64])
+def test_variable_init_relative(fieldset, mode, coord_type, npart=10):
     """Test that checks relative initialisation of custom variables"""
+    lonlat_type = np.float64 if coord_type == 'double' else np.float32
+
     class TestParticle(ptype[mode]):
-        p_base = Variable('p_base', dtype=np.float32, initial=10.)
-        p_relative = Variable('p_relative', dtype=np.float32,
+        p_base = Variable('p_base', dtype=lonlat_type, initial=10.)
+        p_relative = Variable('p_relative', dtype=lonlat_type,
                               initial=attrgetter('p_base'))
-        p_offset = Variable('p_offset', dtype=np.float32,
+        p_offset = Variable('p_offset', dtype=lonlat_type,
                             initial=attrgetter('p_base'))
-        p_lon = Variable('p_lon', dtype=np.float32,
+        p_lon = Variable('p_lon', dtype=lonlat_type,
                          initial=attrgetter('lon'))
-        p_lat = Variable('p_lat', dtype=np.float32,
+        p_lat = Variable('p_lat', dtype=lonlat_type,
                          initial=attrgetter('lat'))
 
         def __init__(self, *args, **kwargs):
             super(TestParticle, self).__init__(*args, **kwargs)
             self.p_offset += 2.
-    lon = np.linspace(0, 1, npart, dtype=np.float32)
-    lat = np.linspace(1, 0, npart, dtype=np.float32)
-    pset = ParticleSet(fieldset, pclass=TestParticle, lon=lon, lat=lat)
+    lon = np.linspace(0, 1, npart, dtype=lonlat_type)
+    lat = np.linspace(1, 0, npart, dtype=lonlat_type)
+    pset = ParticleSet(fieldset, pclass=TestParticle, lon=lon, lat=lat, lonlatdepth_dtype=coord_type)
     # Adjust base variable to test for aliasing effects
     for p in pset:
         p.p_base += 3.
diff --git a/tests/test_scripts.py b/tests/test_scripts.py
index 31fc87d99c..7a473b4954 100644
--- a/tests/test_scripts.py
+++ b/tests/test_scripts.py
@@ -16,7 +16,7 @@ def create_outputfiles(dir):
     endtime = delta(hours=24)
     x = 3. * (1. / 1.852 / 60)
     y = (fieldset.U.lat[0] + x, fieldset.U.lat[-1] - x)
-    lat = np.linspace(y[0], y[1], npart, dtype=np.float32)
+    lat = np.linspace(y[0], y[1], npart)
 
     fp = dir.join("DelayParticle.nc")
     output_file = pset.ParticleFile(name=fp, outputdt=delaytime)