diff --git a/examples/example_FTLE_norkyst.py b/examples/example_FTLE_norkyst.py
new file mode 100644
index 000000000..b58f66941
--- /dev/null
+++ b/examples/example_FTLE_norkyst.py
@@ -0,0 +1,43 @@
+#!/usr/bin/env python
+"""
+LCS Norkyst
+==================================
+"""
+
+from datetime import datetime, timedelta
+import numpy as np
+from opendrift.readers import reader_netCDF_CF_generic
+from opendrift.models.oceandrift import OceanDrift
+
+o = OceanDrift(loglevel=20)  # Set loglevel to 0 for debug information
+
+# Norkyst ocean model
+reader_norkyst = reader_netCDF_CF_generic.Reader(o.test_data_folder() + '16Nov2015_NorKyst_z_surface/norkyst800_subset_16Nov2015.nc')
+
+#o.add_reader([reader_norkyst, reader_arome])
+o.add_reader([reader_norkyst])
+
+
+#%%
+# Calculating attracting/backwards FTLE at 20 hours
+lcs = o.calculate_ftle(
+    time=reader_norkyst.start_time + timedelta(hours=24),
+    time_step=timedelta(minutes=15),
+    duration=timedelta(hours=3), delta=800,
+    RLCS=False)
+
+#%%
+# Simulation from beginning and up to 30 hours (time of LCS)
+o.reset()
+lons = np.linspace(3.2, 5.0, 100)
+lats = np.linspace(59.8, 61, 100)
+lons, lats = np.meshgrid(lons, lats)
+lons = lons.ravel()
+lats = lats.ravel()
+o.seed_elements(lons, lats, radius=0, number=10000,
+                time=reader_norkyst.start_time)
+
+o.run(end_time=reader_norkyst.start_time+timedelta(hours=24),
+      time_step=timedelta(minutes=30))
+
+o.plot(lcs=lcs, vmin=1e-7, vmax=1e-4, cmap='Reds', markersize=1, colorbar=True, show_particles=True, show_initial=False, linewidth=0)
diff --git a/examples/example_LCS_norkyst.py b/examples/example_LCS_norkyst.py
index b3eea0e4f..5ab669f34 100755
--- a/examples/example_LCS_norkyst.py
+++ b/examples/example_LCS_norkyst.py
@@ -5,9 +5,11 @@
 """
 
 from datetime import datetime, timedelta
-
 from opendrift.readers import reader_netCDF_CF_generic
 from opendrift.models.oceandrift import OceanDrift
+import matplotlib.pyplot as plt
+import numpy as np
+import cartopy.crs as ccrs
 
 o = OceanDrift(loglevel=20)  # Set loglevel to 0 for debug information
 
@@ -19,20 +21,109 @@
 
 
 #%%
-# Calculating attracting/backwards FTLE/LCS at 20 hours
-lcs = o.calculate_ftle(
-    time=reader_norkyst.start_time + timedelta(hours=20),
-    time_step=timedelta(minutes=30),
-    duration=timedelta(hours=5), delta=1000,
-    RLCS=False)
+
+x0,y0 = reader_norkyst.lonlat2xy(4.5,60.5)
+d = 30000
+lcsproj = reader_norkyst.proj
+lcs = o.calculate_lcs(
+    time=reader_norkyst.start_time + timedelta(hours=24),
+    time_step=timedelta(minutes=15), reader=lcsproj,
+    duration=timedelta(hours=3), delta=400, domain=[x0-d,x0+d,y0-d,y0+d])
+
+l1 = lcs['eigval'][0,:,:,0]
+l2 = lcs['eigval'][0,:,:,1]
+mask = l2>l1
+
+lmax = l1.copy()
+lmax[mask] = l2[mask]
+
+lmin = l2.copy()
+lmin[mask] = l1[mask]
+
+xi1 = lcs['eigvec'][0,:,:,0]
+xi2 = lcs['eigvec'][0,:,:,1]
+
+stretch = xi1.copy()
+stretch[mask] = xi2[mask]
+
+shrink = xi2.copy()
+shrink[mask] = xi1[mask]
+
+
+# find local maxima of largest eigenvalue
+from skimage.feature import peak_local_max
+#peaks = peak_local_max(lcs['eigval'][0,:,:,1],5)
+peaks = peak_local_max(-lmin,5)
+plon = [lcs['lon'][peaks[i,0],peaks[i,1]] for i in range(peaks.shape[0])]
+plat = [lcs['lat'][peaks[i,0],peaks[i,1]] for i in range(peaks.shape[0])]
+
+
+from opendrift.readers.reader_constant_2d import Reader
+x,y = reader_norkyst.lonlat2xy(lcs['lon'],lcs['lat'])
+r = Reader(x=x, y=y, proj4=reader_norkyst.proj4, array_dict = {'x_sea_water_velocity': stretch[:,:,0], 'y_sea_water_velocity': stretch[:,:,1]})
+
+
+xi = OceanDrift(loglevel=20)
+xi.add_reader(r)
+xi.seed_elements(plon, plat, time=datetime.now())
+xi.run(duration=timedelta(hours=3), time_step=300)
+str_lon = xi.history['lon']
+str_lat = xi.history['lat']
+
+#xi.plot(linewidth=2, linecolor='r',show_particles=False)
+
+xi.reset()
+xi.seed_elements(plon, plat, time=datetime.now())
+xi.run(duration=timedelta(hours=3), time_step=-300)
+str_lon = np.concatenate((str_lon, xi.history['lon'][:,::-1]), axis=0)
+str_lat = np.concatenate((str_lat, xi.history['lat'][:,::-1]), axis=0)
+
 
 #%%
 # Simulation from beginning and up to 30 hours (time of LCS)
+
 o.reset()
-o.seed_elements(lon=4.4, lat=60.2, number=1000, radius=1000,
+lons = np.linspace(4.0, 5.0, 100)
+lats = np.linspace(60., 61., 100)
+lons, lats = np.meshgrid(lons, lats)
+lons = lons.ravel()
+lats = lats.ravel()
+o.seed_elements(lons, lats, radius=0, number=10000,
                 time=reader_norkyst.start_time)
-o.run(end_time=reader_norkyst.start_time+timedelta(hours=20),
+
+o.run(end_time=reader_norkyst.start_time+timedelta(hours=24),
       time_step=timedelta(minutes=30))
+ax, plt = o.plot(cmap='Reds', vmax=2, markersize=1, colorbar=True, show_particles=True, show_initial=False, linewidth=0, show=False)
+
+'''
+fig = plt.figure()
+#proj=reader_norkyst.proj
+
+lon, lat = lcs['lon'], lcs['lat']
+x,y = reader_norkyst.lonlat2xy(lon, lat)
+px,py = reader_norkyst.lonlat2xy(plon, plat)
+#strx, stry = reader_norkyst.lonlat2xy(str_lon, str_lat)
+fig.add_subplot(121)
+plt.pcolormesh(x, y, np.log(np.sqrt(lmin)),vmin=-1,vmax=3, cmap='cividis'),plt.colorbar()
+plt.quiver(x[::5,::5], y[::5,::5], shrink[::5,::5,0], shrink[::5,::5,1])
+plt.plot(px, py, '.r')
+plt.title('shrink lines')
+fig.add_subplot(122)
+
+plt.pcolormesh(x,y , np.log(np.sqrt(lmax)), cmap='cividis'),plt.colorbar()
+plt.quiver(x[::5,::5], y[::5,::5], stretch[::5,::5,0], stretch[::5,::5,1])
+
+plt.title('stretch lines')
+
+<<<<<<< HEAD
+plt.figure()
+'''
+
+gcrs = ccrs.PlateCarree()
+ax.plot(str_lon.T, str_lat.T, '-r', ms=0.5, transform=gcrs)
+plt.show()
 
 o.plot(lcs=lcs, vmin=1e-7, vmax=1e-4, colorbar=True, show_elements=True)
 
+
+#o.animation(buffer=0, lcs=ftle, hide_landmask=True)
diff --git a/examples/example_double_gyre_LCS_particles.py b/examples/example_double_gyre_LCS_particles.py
index b9483ae18..95e1f7ea7 100755
--- a/examples/example_double_gyre_LCS_particles.py
+++ b/examples/example_double_gyre_LCS_particles.py
@@ -3,7 +3,7 @@
 Double gyre - LCS with particles
 ============================================
 
-Drift of particles in an idealised (analytical) eddy current field, 
+Drift of particles in an idealised (analytical) eddy current field,
 plotted on top of the LCS. This takes some minutes to calculate.
 """
 
@@ -53,8 +53,7 @@
 o.disable_vertical_motion()
 o.run(duration=duration, time_step=time_step,
       time_step_output=time_step_output)
-o.animation(buffer=0, lcs=lcs, hide_landmask=True)
+o.animation(buffer=0, lcs=lcs, cmap='cividis', hide_landmask=True)
 
 #%%
 # .. image:: /gallery/animations/example_double_gyre_LCS_particles_0.gif
-
diff --git a/examples/example_double_gyre_LCS_snapshot.py b/examples/example_double_gyre_LCS_snapshot.py
new file mode 100644
index 000000000..59039611a
--- /dev/null
+++ b/examples/example_double_gyre_LCS_snapshot.py
@@ -0,0 +1,75 @@
+#!/usr/bin/env python
+"""
+Double gyre - LCS with particles
+============================================
+
+Drift of particles in an idealised (analytical) eddy current field,
+plotted on top of the LCS. This takes some minutes to calculate.
+"""
+
+from datetime import datetime, timedelta
+import matplotlib.pyplot as plt
+import numpy as np
+
+from opendrift.readers import reader_double_gyre
+from opendrift.models.oceandrift import OceanDrift
+
+#%%
+# Setting some parameters
+plot_time =  timedelta(seconds=12)
+duration = timedelta(seconds=12)  # T
+time_step=timedelta(seconds=.5)
+time_step_output=timedelta(seconds=.5)
+delta=.01  # spatial resolution
+#steps = int(duration.total_seconds()/ time_step_output.total_seconds() + 1)
+
+o = OceanDrift(loglevel=20)
+
+#%%
+# Note that Runge-Kutta here makes a difference to Euler scheme
+o.set_config('drift:scheme', 'runge-kutta4')
+o.disable_vertical_motion()
+o.fallback_values['land_binary_mask'] = 0
+
+double_gyre = reader_double_gyre.Reader(epsilon=.25, omega=0.628, A=0.1)
+print(double_gyre)
+o.add_reader(double_gyre)
+
+
+#%%
+# Calculate Lyapunov exponents
+#times = [double_gyre.initial_time +  n*time_step_output for n in range(steps)]
+ftle = o.calculate_ftle(time=double_gyre.initial_time + plot_time, time_step=time_step,
+                       duration=duration, delta=delta, RLCS=False)
+
+lcs = o.calculate_lcs(time=double_gyre.initial_time + plot_time, time_step=-time_step,
+                       duration=duration, delta=delta)
+#%%
+# Make run with particles for the same period
+o.reset()
+x = [.9]
+y = [.5]
+lon, lat = double_gyre.xy2lonlat(x, y)
+
+o.seed_elements(lon, lat, radius=.15, number=2000,
+                time=double_gyre.initial_time)
+
+o.disable_vertical_motion()
+o.run(duration=plot_time, time_step=time_step,
+      time_step_output=time_step_output)
+lonmin, latmin = double_gyre.xy2lonlat(0.,0.)
+lonmax, latmax = double_gyre.xy2lonlat(2.,1.)
+o.plot(lcs=ftle, show_initial=False, linewidth = 0, corners =[lonmin, lonmax, latmin, latmax], cmap='cividis')
+
+fig = plt.figure()
+fig.add_subplot(211)
+plt.pcolormesh(np.log(np.sqrt(lcs['eigval'][0,:,:,0])), cmap='cividis'),plt.colorbar()
+plt.quiver(lcs['eigvec'][0,:,:,0,0], lcs['eigvec'][0,:,:,0,1])
+fig.add_subplot(212)
+plt.pcolormesh(np.log(np.sqrt(lcs['eigval'][0,:,:,1])), cmap='cividis'),plt.colorbar()
+plt.quiver(lcs['eigvec'][0,:,:,0,0], lcs['eigvec'][0,:,:,0,1])
+plt.title('Eigenvalues and Eigenvectors of Cauchy-Green Strain tensor')
+#o.animation(buffer=0, lcs=ftle, hide_landmask=True)
+
+#%%
+# .. image:: /gallery/animations/example_double_gyre_LCS_particles_0.gif
diff --git a/opendrift/models/basemodel.py b/opendrift/models/basemodel.py
index 52c58415c..4a36fd93b 100644
--- a/opendrift/models/basemodel.py
+++ b/opendrift/models/basemodel.py
@@ -1829,11 +1829,11 @@ def seed_within_polygon(self, lons, lats, number=None, **kwargs):
         area = 0.0
         for i in range(-1, len(x)-1):
             area += x[i] * (y[i+1] - y[i-1])
-        area = abs(area) / 2
 
         # Make points, evenly distributed
         deltax = np.sqrt(area/number)
         lonpoints = np.array([])
+        area = abs(area) / 2
         latpoints = np.array([])
         lonlat = poly.get_xy()
         lon = lonlat[:, 0]
@@ -2131,7 +2131,7 @@ def set_fallback_values(self, refresh=False):
 
     def run(self, time_step=None, steps=None, time_step_output=None,
             duration=None, end_time=None, outfile=None, export_variables=None,
-            export_buffer_length=100, stop_on_error=False):
+            export_buffer_length=100, stop_on_error=False, move_from_land=True):
         """Start a trajectory simulation, after initial configuration.
 
         Performs the main loop:
@@ -2357,15 +2357,17 @@ def run(self, time_step=None, steps=None, time_step_output=None,
 
             self.timer_end('preparing main loop:making dynamical landmask')
 
-        # Move point seeded on land to ocean
-        if self.get_config('seed:ocean_only') is True and \
-            ('land_binary_mask' in self.required_variables):
-            #('land_binary_mask' not in self.fallback_values) and \
-            self.timer_start('preparing main loop:moving elements to ocean')
-            self.elements_scheduled.lon, self.elements_scheduled.lat = \
+
+        # Move point seed on land to ocean
+        if move_from_land==True:
+            if self.get_config('seed:ocean_only') is True and \
+                ('land_binary_mask' not in self.fallback_values) and \
+                ('land_binary_mask' in self.required_variables):
+                self.timer_start('preparing main loop:moving elements to ocean')
+                self.elements_scheduled.lon, self.elements_scheduled.lat = \
                 self.closest_ocean_points(self.elements_scheduled.lon,
                                           self.elements_scheduled.lat)
-            self.timer_end('preparing main loop:moving elements to ocean')
+                self.timer_end('preparing main loop:moving elements to ocean')
 
         ####################################################################
         # Preparing history array for storage in memory and eventually file
@@ -4611,8 +4613,10 @@ def calculate_ftle(self, reader=None, delta=None, domain=None,
             # Backwards
             if ALCS is True:
                 self.reset()
+                #self.seed_elements(lons.ravel(), lats.ravel(),
+                #                   time=t+duration, z=z)
                 self.seed_elements(lons.ravel(), lats.ravel(),
-                                   time=t+duration, z=z)
+                                   time=t, z=z)
                 self.run(duration=duration, time_step=-time_step)
                 b_x1, b_y1 = proj(
                     self.history['lon'].T[-1].reshape(X.shape),
@@ -4626,6 +4630,71 @@ def calculate_ftle(self, reader=None, delta=None, domain=None,
 
         return lcs
 
+    def calculate_lcs(self, reader=None, delta=None, domain=None,
+                       time=None, time_step=None, duration=None, z=0):
+
+        if reader is None:
+            self.logger.info('No reader provided, using first available:')
+            reader = list(self.readers.items())[0][1]
+            self.logger.info(reader.name)
+        if isinstance(reader, pyproj.Proj):
+            proj = reader
+        elif isinstance(reader,str):
+            proj = pyproj.Proj(reader)
+        else:
+            proj = reader.proj
+
+        import scipy.ndimage as ndimage
+        from opendrift.models.physics_methods import cg_eigenvectors
+
+        if not isinstance(duration, timedelta):
+            duration = timedelta(seconds=duration)
+
+        if domain==None:
+            xs = np.arange(reader.xmin, reader.xmax, delta)
+            ys = np.arange(reader.ymin, reader.ymax, delta)
+        else:
+            xmin, xmax, ymin, ymax = domain
+            xs = np.arange(xmin, xmax, delta)
+            ys = np.arange(ymin, ymax, delta)
+
+        X, Y = np.meshgrid(xs, ys)
+        lons, lats = proj(X, Y, inverse=True)
+
+        if time is None:
+            time = reader.start_time
+        if not isinstance(time, list):
+            time = [time]
+        # dictionary to hold LCS calculation
+        lcs = {'time': time, 'lon': lons, 'lat':lats}
+        lcs['LCS'] = np.zeros((len(time), len(ys), len(xs)))
+        lcs['eigvec'] = np.zeros((len(time), len(ys), len(xs), 2, 2 ))
+        lcs['eigval'] = np.zeros((len(time), len(ys), len(xs), 2))
+
+        T = np.abs(duration.total_seconds())
+        for i, t in enumerate(time):
+            self.logger.info('Calculating LCS for ' + str(t))
+            # Forward or bachward flow map
+            self.reset()
+            self.seed_elements(lons.ravel(), lats.ravel(),
+                               time=t, z=z)
+            self.run(duration=duration, time_step=time_step, move_from_land=False)
+            x1, y1 = proj(
+                self.history['lon'].T[-1].reshape(X.shape),
+                self.history['lat'].T[-1].reshape(X.shape))
+            lamba,xi =  cg_eigenvectors(x1-X, y1-Y, delta, T)
+            lcs['eigval'][i,:,:,:] = lamba
+            lcs['eigvec'][i,:,:,:,:] = xi
+
+        lcs['eigvec'] = np.ma.masked_invalid(lcs['eigvec'])
+        lcs['eigval'] = np.ma.masked_invalid(lcs['eigval'])
+        # Flipping ALCS left-right if using backward time flow map; Not sure why this is needed
+        if time_step.total_seconds() < 0:
+            lcs['eigval'] = lcs['eigval'][:,::-1,::-1]
+            lcs['eigvec'] = lcs['eigvec'][:,::-1,::-1]
+        return lcs
+
+
     def center_of_gravity(self, onlysurface=False):
         """
         calculate center of mass and variance of all elements
diff --git a/opendrift/models/lophelia.py b/opendrift/models/lophelia.py
new file mode 100644
index 000000000..c5d7a5c8a
--- /dev/null
+++ b/opendrift/models/lophelia.py
@@ -0,0 +1,180 @@
+# This file is part of OpenDrift.
+#
+# OpenDrift is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, version 2
+#
+# OpenDrift is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with OpenDrift.  If not, see <https://www.gnu.org/licenses/>.
+#
+# Copyright 2015, Knut-Frode Dagestad, MET Norway
+
+import numpy as np
+import logging; logger = logging.getLogger(__name__)
+
+from opendrift.models.oceandrift import OceanDrift, Lagrangian3DArray
+#from opendrift.elements import LagrangianArray
+
+
+# Defining element properties
+class LopheliaLarvae(Lagrangian3DArray):
+    """Extending Lagrangian3DArray with specific properties for Lophelia larvae
+    """
+
+    variables = Lagrangian3DArray.add_variables([
+        ('diameter', {'dtype': np.float32,
+                      'units': 'm',
+                      'default': 0.0002}),  # For Lophelia
+        ('density', {'dtype': np.float32,
+                     'units': 'kg/m^3',
+                     'default': 1028.}),
+        ('age_seconds', {'dtype': np.float32,
+                         'units': 's',
+                         'default': 0.}),
+        ('competence', {'dtype': np.float32,
+                         'units': 's',
+                         'default': 0.}),
+        ('hatched', {'dtype': np.float32,
+                     'units': '',
+                     'default': 0.})])
+
+class LopheliaLarvaeDrift(OceanDrift):
+    """Buoyant particle trajectory model based on the OpenDrift framework.
+
+        Developed at MET Norway
+
+        Generic module for particles that are subject to vertical turbulent
+        mixing with the possibility for positive or negative buoyancy
+
+        Lophelia larvae assumed to ber neutrally buoyant throughout
+
+        Particles could be e.g. oil droplets, plankton, or sediments
+
+        Under construction.
+    """
+
+
+
+    ElementType = LopheliaLarvae
+
+    required_variables = {
+        'x_sea_water_velocity': {'fallback': 0},
+        'y_sea_water_velocity': {'fallback': 0},
+        'sea_surface_wave_stokes_drift_x_velocity': {'fallback': 0},
+        'sea_surface_wave_stokes_drift_y_velocity': {'fallback': 0},
+        'sea_surface_wave_significant_height': {'fallback': 0},
+        'sea_ice_area_fraction': {'fallback': 0},
+        'x_wind': {'fallback': 0},
+        'y_wind': {'fallback': 0},
+        'land_binary_mask': {'fallback': None},
+        'sea_floor_depth_below_sea_level': {'fallback': 100},
+        'ocean_vertical_diffusivity': {'fallback': 0.02, 'profiles': True},
+        'ocean_mixed_layer_thickness': {'fallback': 50},
+        'sea_water_temperature': {'fallback': 10, 'profiles': True},
+        'sea_water_salinity': {'fallback': 34, 'profiles': True},
+        'surface_downward_x_stress': {'fallback': 0},
+        'surface_downward_y_stress': {'fallback': 0},
+        'turbulent_kinetic_energy': {'fallback': 0},
+        'turbulent_generic_length_scale': {'fallback': 0},
+        'upward_sea_water_velocity': {'fallback': 0},
+      }
+
+    # Vertical profiles of the following parameters will be available in
+    # dictionary self.environment.vertical_profiles
+    # E.g. self.environment_profiles['x_sea_water_velocity']
+    # will be an array of size [vertical_levels, num_elements]
+    # The vertical levels are available as
+    # self.environment_profiles['z'] or
+    # self.environment_profiles['sigma'] (not yet implemented)
+    required_profiles = ['sea_water_temperature',
+                         'sea_water_salinity',
+                         'ocean_vertical_diffusivity']
+    # The depth range (in m) which profiles shall cover
+    required_profiles_z_range = [-3000, 0]
+
+
+    # Default colors for plotting
+    status_colors = {'initial': 'green', 'active': 'blue',
+                     'spawned': 'red', 'settled': 'yellow', 'died': 'magenta'}
+
+
+    def __init__(self, *args, **kwargs):
+
+        # Calling general constructor of parent class
+        super(LopheliaLarvaeDrift, self).__init__(*args, **kwargs)
+
+        # By default, elements do not strand towards coastline
+        self.set_config('general:coastline_action', 'previous')
+
+        # Vertical mixing is enabled by default
+        self.set_config('drift:vertical_mixing', True)
+
+    def update_terminal_velocity(self, z_index=None):
+        """Calculate terminal velocity for larvae
+
+        """
+        g = 9.81  # ms-2
+
+        # Properties that determine element buoyancy
+        larvaesize = self.elements.diameter  # 0.0002 for Lophelia
+        #eggsalinity = self.elements.neutral_buoyancy_salinity
+
+
+	    # Set vertical velocities depending on stage (age) in larval phase:
+
+        scenario = 'short' # Set scenario, 'short' or 'long' depending on length of pre-competency period
+        sday = 24*60*60 # Seconds in one day
+        self.reftemp = 8 # degree C
+
+        if scenario == 'short':
+            dup2dwn = 21
+        elif scenario == 'long':
+    	    dup2dwn = 42
+
+        idd = np.arange(15, dtype=float)
+        idd = np.insert(idd, [7, 8], [6.001, 7.001])
+        idd = np.append(idd, [dup2dwn, dup2dwn+1])
+        idd = idd*sday*self.reftemp
+
+        w_idd = np.array([0, 0, 0.00005, 0.00005, 0.00005, 0.0002, 0.00025, 0, 0, 0, 0.0003, 0.00035, 0.0004, 0.00045, 0.0005, 0.00055, 0.0006, 0.0006, -0.0006]) # Vertical velocities at specified stages
+
+        W = np.interp(self.elements.competence, idd, w_idd)
+
+        print(w)
+        self.elements.terminal_velocity = W
+
+    def update(self):
+        """Update positions and properties of particles."""
+
+        # Stokes drift
+        self.stokes_drift()
+
+        # Update element age
+        self.elements.age_seconds += self.time_step.total_seconds()
+
+        # Update competence
+        self.elements.competence += self.time_step.total_seconds() * self.environment.sea_water_temperature
+
+        # Turbulent Mixing
+        self.update_terminal_velocity()
+        self.vertical_mixing()
+        #self.vertical_buoyancy()
+
+        # Horizontal advection
+        self.advect_ocean_current()
+
+        # Vertical advection
+        if self.get_config('drift:vertical_advection') is True:
+            self.vertical_advection()
+
+        if self.time_step.total_seconds() < 0:
+            self.deactivate_elements(self.elements.competence <= 0., reason='spawned')
+
+        if self.time_step.total_seconds() > 0:
+            self.deactivate_elements(self.elements.competence > 21.*24*3600* self.reftemp & self.elements.z < -sea_floor_depth, reason='settled')
+            self.deactivate_elements(self.elements.competence > 60.*24*3600* self.reftemp , reason='died')
diff --git a/opendrift/models/physics_methods.py b/opendrift/models/physics_methods.py
index 54551633c..1d2d75b69 100644
--- a/opendrift/models/physics_methods.py
+++ b/opendrift/models/physics_methods.py
@@ -277,7 +277,7 @@ def stokes_drift_profile_breivik(stokes_u_surface, stokes_v_surface,
 
 
 def ftle(X, Y, delta, duration):
-    """Calculate Finite Time Lyapunov Exponents"""
+    """Calculate Finite Time Lyapunov Exponents from flow map"""
     # From Johannes Rohrs
     nx = X.shape[0]
     ny = X.shape[1]
@@ -304,6 +304,38 @@ def ftle(X, Y, delta, duration):
 
     return FTLE
 
+def cg_eigenvectors(X, Y, delta, duration):
+    """Calculate eigenvector and eigenvalues of cauchy-green strain tensor
+    returns eigenvalues, eigenvectors of the flow map X,Y
+    """
+
+    # From Johannes Rohrs
+    nx = X.shape[0]
+    ny = X.shape[1]
+    J = np.empty([nx,ny,2,2],np.float)
+    lamba = np.empty([nx,ny,2],np.float)
+    xi = np.empty([nx,ny,2,2],np.float)
+
+    # gradient
+    dx = np.gradient(X)
+    dy = np.gradient(Y)
+
+    # Jacobian
+    J[:,:,0,0] = dx[0] / (2*delta)
+    J[:,:,1,0] = dy[0] / (2*delta)
+    J[:,:,0,1] = dx[1] / (2*delta)
+    J[:,:,1,1] = dy[1] / (2*delta)
+
+    for i in range(0,nx):
+        for j in range(0,ny):
+            # Green-Cauchy tensor
+            D = np.dot(np.transpose(J[i,j]), J[i,j])
+            # its largest eigenvalue
+            lamba[i,j], xi[i,j] = np.linalg.eig(D)
+            lamba[lamba>1000]=np.nan # mask out land artifacts
+    return lamba, xi
+
+
 class PhysicsMethods:
     """Physics methods to be inherited by OpenDriftSimulation class"""