diff --git a/aeolis/constants.py b/aeolis/constants.py
index 5c50e094..137f64be 100644
--- a/aeolis/constants.py
+++ b/aeolis/constants.py
@@ -322,6 +322,7 @@
     'method_transport'              : 'bagnold',          # Name of method to compute equilibrium sediment transport rate
     'method_roughness'              : 'constant',         # Name of method to compute the roughness height z0, note that here the z0 = k, which does not follow the definition of Nikuradse where z0 = k/30.
     'method_grainspeed'             : 'windspeed',        # Name of method to assume/compute grainspeed (windspeed, duran, constant)
+    'method_shear'                  : 'fft',              # Name of method to compute topographic effects on wind shear stress (fft, quasi2d, duna2d (experimental))
     'max_error'                     : 1e-6,               # [-] Maximum error at which to quit iterative solution in implicit numerical schemes
     'max_iter'                      : 1000,               # [-] Maximum number of iterations at which to quit iterative solution in implicit numerical schemes
     'max_iter_ava'                  : 1000,               # [-] Maximum number of iterations at which to quit iterative solution in avalanching calculation
diff --git a/aeolis/rotation.py b/aeolis/rotation.py
new file mode 100644
index 00000000..fcce4b23
--- /dev/null
+++ b/aeolis/rotation.py
@@ -0,0 +1,306 @@
+'''This file is part of AeoLiS.
+
+AeoLiS 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, either version 3 of the License, or
+(at your option) any later version.
+
+AeoLiS 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 AeoLiS.  If not, see <http://www.gnu.org/licenses/>.
+
+AeoLiS  Copyright (C) 2015 Bas Hoonhout
+
+bas.hoonhout@deltares.nl         b.m.hoonhout@tudelft.nl
+Deltares                         Delft University of Technology
+Unit of Hydraulic Engineering    Faculty of Civil Engineering and Geosciences
+Boussinesqweg 1                  Stevinweg 1
+2629 HVDelft                     2628CN Delft
+The Netherlands                  The Netherlands
+
+'''
+
+from __future__ import absolute_import, division
+import logging
+import aeolis.vegetation
+from aeolis.utils import *
+
+# initialize logger
+logger = logging.getLogger(__name__)
+
+
+class rotationClass:
+    '''XXXX.
+    '''
+
+    igrid = {}
+    cgrid = {}
+    istransect = False
+
+    def __init__(self, x, y, z, dx, dy, buffer_width):
+        '''Class initialization
+
+        Parameters
+        ----------
+        x : numpy.ndarray
+            2D array with x-coordinates of input grid
+        y : numpy.ndarray
+            2D array with y-coordinates of input grid
+        z : numpy.ndarray
+            2D array with topography of input grid
+        dx : float, optional
+            Grid spacing in x dimension of computational grid
+            (default: 1)
+        dy : float, optional
+            Grid spacing of y dimension of computational grid
+            (default: 1)
+        buffer_width : float, optional
+            Width of buffer distance between input grid boundary and
+            computational grid boundary (default: 100)
+        '''
+
+        if z.shape[0] == 1:
+            self.istransect = True
+
+        # Assigning values to original (i) and computational (c) grid
+        self.cgrid = dict(dx=dx, dy=dy)
+
+        # Setting buffer settings
+        self.buffer_width = buffer_width
+
+    def __call__(self, x, y, z, udir, ustar, tau, hveg, type_flag, params):
+        # Reload x and y because of horizontalized input-grid
+        self.igrid = dict(x=x, y=y, z=z, ustar=ustar, hveg=hveg, tau=tau)
+
+        # Convert to cartesian to perform all the rotations
+        u_angle = 270. - udir  # wind angle
+
+        # Creating the computational grid
+        self.set_computational_grid(udir)
+
+        # Storing computational (c) and original (i) grids
+        gi = self.igrid  # initial grid
+        gc = self.cgrid  # computational grid
+
+        # Rotate computational (c) grid to the current wind direction
+        gc['x'], gc['y'] = self.rotate(gc['xi'], gc['yi'], -u_angle, origin=(self.x0, self.y0))
+        # Interpolate bed levels and shear to the computational grid
+        gc['z'] = self.interpolate(gi['x'], gi['y'], gi['z'], gc['x'], gc['y'], 0)
+        gc['x'], gc['y'] = self.rotate(gc['x'], gc['y'], u_angle, origin=(self.x0, self.y0))
+        gi['x'], gi['y'] = self.rotate(gi['x'], gi['y'], u_angle, origin=(self.x0, self.y0))
+
+        #Do relevant calculations here
+        if type_flag == 1: #okin calcs
+            gc['ustar'] = self.interpolate(gi['x'], gi['y'], gi['ustar'], gc['x'], gc['y'], 0)
+            gc['hveg'] = self.interpolate(gi['x'], gi['y'], gi['hveg'], gc['x'], gc['y'], 0)
+            gc['ustar'] = aeolis.vegetation.compute_okin_shear(gc['x'], gc['ustar'], gc['hveg'], params)
+            gi['ustar'] = self.interpolate(gc['x'], gc['y'], gc['ustar'], gi['x'], gi['y'], 0)
+        if type_flag == 2: #sandfence calcs
+            gc['ustar'] = self.interpolate(gi['x'], gi['y'], gi['ustar'], gc['x'], gc['y'], 0)
+            gc['fence_height'] = self.interpolate(gi['x'], gi['y'], gi['fence_height'], gc['x'], gc['y'], 0)
+            gc['ustar'] = aeolis.fences.compute_fence_shear(gc['x'], gc['ustar'], gc['fence_height'], params)
+            gi['ustar'] = self.interpolate(gc['x'], gc['y'], gc['ustar'], gi['x'], gi['y'], 0)
+        if type_flag == 3: #shear perturbation calcs
+            gc['tau'] = self.interpolate(gi['x'], gi['y'], gi['tau'], gc['x'], gc['y'], 0)
+            if params['shear_type'] == 'duna2d':
+                gc['tau'] = aeolis.wind.compute_shear_perturbation_1D(gc['x'], gc['z'], gc['tau'], params)
+            if params['shear_type'] == 'quasi2d':
+                gc['tau'] = aeolis.wind.compute_shear_perturbation_kroy1D(gc['x'], gc['y'], gc['z'], gc['tau'], params)
+            gc['zsep'], gc['hsep'], gc['tau'] = aeolis.wind.separation_quasi2d(gc['x'], gc['y'], gc['z'], gc['tau'], gc['dx'], gc['dy'], params)
+            gi['zsep'] = self.interpolate(gc['x'], gc['y'], gc['zsep'], gi['x'], gi['y'], 0)
+            gi['hsep'] = self.interpolate(gc['x'], gc['y'], gc['hsep'], gi['x'], gi['y'], 0)
+            gi['tau'] = self.interpolate(gc['x'], gc['y'], gc['tau'], gi['x'], gi['y'], 0)
+
+        self.igrid = gi # initial grid
+        self.cgrid = gc  # computational grid
+
+        return self
+
+    # Input functions for __call()
+    def set_computational_grid(self, udir):
+        '''Define computational grid
+
+        The computational grid is square with dimensions equal to the
+        diagonal of the bounding box of the input grid, plus twice the
+        buffer width.
+
+        '''
+
+        # Copying the original (i) and computational (c) grid
+        gi = self.igrid
+        gc = self.cgrid
+
+        # Compute grid center, same for both original (i) and computational (c) grid
+        x0, y0 = np.mean(gi['x']), np.mean(gi['y'])
+
+        # Initialization
+        b_W = np.zeros(4)
+        b_L = np.zeros(4)
+        xcorner = np.zeros(4)
+        ycorner = np.zeros(4)
+
+        # Computing the corner-points of the grid
+        xcorner[0] = gi['x'][0, 0]
+        ycorner[0] = gi['y'][0, 0]
+        xcorner[1] = gi['x'][-1, 0]
+        ycorner[1] = gi['y'][-1, 0]
+        xcorner[2] = gi['x'][0, -1]
+        ycorner[2] = gi['y'][0, -1]
+        xcorner[3] = gi['x'][-1, -1]
+        ycorner[3] = gi['y'][-1, -1]
+
+        # Preventing vertical lines
+        udir_verticals = np.arange(-1080, 1080, 90)
+        udir_vertical_bool = False
+        for udir_vertical in udir_verticals:
+            if (abs(udir - udir_vertical) <= 0.001):
+                udir_vertical_bool = True
+        if udir_vertical_bool:
+            udir -= 0.1
+
+        # Compute slope (m) and intercept (b) from parallel lines along all (4) grids corners
+        for i in range(4):
+            # Parallel boundaries
+            m_W, b_W[i] = np.polyfit([xcorner[i], xcorner[i] - np.sin(np.deg2rad(udir))],
+                                     [ycorner[i], ycorner[i] - np.cos(np.deg2rad(udir))], 1)
+            # Perpendicular boundaries
+            m_L, b_L[i] = np.polyfit([xcorner[i], xcorner[i] - np.sin(np.deg2rad(udir - 90.))],
+                                     [ycorner[i], ycorner[i] - np.cos(np.deg2rad(udir - 90.))], 1)
+
+        # Determine the most outer boundaries (for parallel and perpendicular)
+        db_W = self.maxDiff(b_W)
+        db_L = self.maxDiff(b_L)
+
+        # Compute the distance between the outer boundaries to determine the width (W) and length (L) of the grid
+        self.Width = abs(db_W) / np.sqrt((m_W ** 2.) + 1) + self.buffer_width * 2.
+        self.Length = abs(db_L) / np.sqrt((m_L ** 2.) + 1) + self.buffer_width * 2.
+
+        # Create the grid
+        xc, yc = self.get_exact_grid(x0 - self.Length / 2., x0 + self.Length / 2.,
+                                     y0 - self.Width / 2., y0 + self.Width / 2.,
+                                     gc['dx'], gc['dy'])
+
+        # Storing grid parameters
+        self.x0 = x0
+        self.y0 = y0
+        gc['xi'] = xc
+        gc['yi'] = yc
+
+        return self
+
+    @staticmethod
+    def get_exact_grid(xmin, xmax, ymin, ymax, dx, dy):
+        '''Returns a grid with given gridsizes approximately within given bounding box'''
+
+        x = np.arange(np.floor(xmin / dx) * dx,
+                      np.ceil(xmax / dx) * dx, dx)
+        y = np.arange(np.floor(ymin / dy) * dy,
+                      np.ceil(ymax / dy) * dy, dy)
+        x, y = np.meshgrid(x, y)
+
+        return x, y
+
+    @staticmethod
+    def get_borders(x):
+        '''Returns borders of a grid as one-dimensional array'''
+
+        return np.concatenate((x[0, :].T,
+                               x[1:-1, -1],
+                               x[-1, ::-1].T,
+                               x[-1:1:-1, 0],
+                               x[0, :1]), axis=0)
+
+    @staticmethod
+    def rotate(x, y, alpha, origin=(0, 0)):
+        '''Rotate a matrix over given angle around given origin'''
+
+        xr = x - origin[0]
+        yr = y - origin[1]
+
+        a = alpha / 180. * np.pi
+
+        R = np.asmatrix([[np.cos(a), -np.sin(a)],
+                         [np.sin(a), np.cos(a)]])
+
+        xy = np.concatenate((xr.reshape((-1, 1)),
+                             yr.reshape((-1, 1))), axis=1) * R
+
+        return (np.asarray(xy[:, 0].reshape(x.shape) + origin[0]),
+                np.asarray(xy[:, 1].reshape(y.shape) + origin[1]))
+
+    def get_ustar(self):
+        '''Returns ustar
+        '''
+
+        ustar = self.igrid['ustar']
+        return ustar
+
+    def get_tau(self):
+        '''Returns ustar
+        '''
+
+        tau = self.igrid['tau']
+        return tau
+
+    def maxDiff(self, arr):
+
+        result = 0
+        n = len(arr)
+
+        # Iterate through all pairs.
+        for i in range(0, n):
+            for j in range(i, n):
+
+                if (abs(arr[i] - arr[j]) + abs(i - j) > result):
+                    result = abs(arr[i] - arr[j]) + abs(i - j)
+
+        return result
+
+    def interpolate(self, x, y, z, xi, yi, z0):
+        '''Interpolate one grid to an other'''
+
+        # First compute angle with horizontal
+        dx = x[0, 1] - x[0, 0]
+        dy = y[0, 1] - y[0, 0]
+
+        angle = np.rad2deg(np.arctan(dy / dx))
+
+        if dx <= 0 and dy <= 0:
+            angle += 180.
+
+        # Rotate grids to allign with horizontal
+        x, y = self.rotate(x, y, angle, origin=(self.x0, self.y0))
+        xi, yi = self.rotate(xi, yi, angle, origin=(self.x0, self.y0))
+
+        # Rotate 180 deg if necessary
+        if not np.all(sorted(y[:, 0]) == y[:, 0]) and not np.all(sorted(x[0, :]) == x[0, :]):
+            x, y = self.rotate(x, y, 180, origin=(self.x0, self.y0))
+            xi, yi = self.rotate(xi, yi, 180, origin=(self.x0, self.y0))
+
+        # Concatenate
+        xy = np.concatenate((y.reshape((-1, 1)),
+                             x.reshape((-1, 1))), axis=1)
+
+        xyi = np.concatenate((yi.reshape((-1, 1)),
+                              xi.reshape((-1, 1))), axis=1)
+
+        # Interpolate
+        pad_w = np.maximum(np.shape(x)[0], np.shape(x)[1])
+        x_pad = np.pad(x, ((pad_w, pad_w), (pad_w, pad_w)), 'reflect', reflect_type='odd')
+        y_pad = np.pad(y, ((pad_w, pad_w), (pad_w, pad_w)), 'reflect', reflect_type='odd')
+        z_pad = np.pad(z, ((pad_w, pad_w), (pad_w, pad_w)), 'edge')
+
+        if self.istransect:
+            zi = np.interp(xi.flatten(), x_pad.flatten(), z_pad.flatten()).reshape(xi.shape)
+        else:
+            # in the scipy 1.10 version the regular grid interpolator does not work with non c-contigous arrays.
+            # Here we make a copy as a dirty solution feeding the interpolator with ordered copies
+            inter = scipy.interpolate.RegularGridInterpolator(
+                (y_pad[:, 0].copy(order='C'), x_pad[0, :].copy(order='C')), z_pad, bounds_error=False, fill_value=z0)
+            zi = inter(xyi).reshape(xi.shape)
+
+        return zi
\ No newline at end of file
diff --git a/aeolis/shear.py b/aeolis/shear.py
index dd33ee52..fb3f7485 100644
--- a/aeolis/shear.py
+++ b/aeolis/shear.py
@@ -877,8 +877,7 @@ def interpolate(self, x, y, z, xi, yi, z0):
             inter = scipy.interpolate.RegularGridInterpolator((y_pad[:,0].copy(order='C'), x_pad[0,:].copy(order='C')), z_pad, bounds_error = False, fill_value = z0)
             zi = inter(xyi).reshape(xi.shape)
             
-            
-            
+        
         return zi
     
 
diff --git a/aeolis/vegetation.py b/aeolis/vegetation.py
index 23ef7044..809c6a86 100644
--- a/aeolis/vegetation.py
+++ b/aeolis/vegetation.py
@@ -31,11 +31,21 @@
 import math
 #import matplotlib.pyplot as plt
 from aeolis.wind import *
+import aeolis.rotation
 
 # package modules
 import aeolis.wind
 #from aeolis.utils import *
 
+import numpy as np
+from scipy.integrate import quad
+from scipy.integrate import romberg
+from scipy.optimize import root_scalar
+from scipy.signal import convolve
+from scipy.special import sici
+from scipy.special import erfi
+import matplotlib.pyplot as plt
+
 # initialize logger
 logger = logging.getLogger(__name__)
 
@@ -76,12 +86,27 @@ def initialize (s,p):
     s['germinate'][:,:] = False
     s['lateral'][:,:] = False
 
-    return s
+    if p['vegshear_type'] == 'okin':
+        s['okin'] = aeolis.rotation.rotationClass(s['x'], s['y'], s['zb'],
+                                        dx=p['dx'], dy=p['dy'],
+                                        buffer_width=100)
 
+    return s
 
 def vegshear(s, p):
-    if p['vegshear_type'] == 'okin' and p['ny'] == 0:
-        s = vegshear_okin(s, p)
+    if p['vegshear_type'] == 'okin':
+
+        okin_params = {
+            'okin_c1_veg': p['okin_c1_veg'],
+            'okin_initialred_veg': p['okin_initialred_veg']
+        }
+
+        s['okin'](x=s['x'], y=s['y'], z=s['zb'], udir=s['udir'][0, 0], ustar=s['ustar'], tau=s['tau'], hveg=s['hveg'], type_flag = 1, params=okin_params)
+
+        s['ustar'] = s['okin'].get_ustar()
+        s['ustars'] = - s['ustar'] * np.sin((-p['alfa'] + s['udir']) / 180. * np.pi)
+        s['ustarn'] = - s['ustar'] * np.cos((-p['alfa'] + s['udir']) / 180. * np.pi)
+
     else:
         s = vegshear_raupach(s, p)
 
@@ -325,3 +350,69 @@ def vegshear_raupach(s, p):
     s['ustarn'] = s['ustar'] * etn
 
     return s
+
+def compute_okin_shear(xgrd, ustar,  hveg, okin_params):
+    #Approach to calculate shear reduction in the lee of plants using the general approach of:
+    #Okin (2008), JGR, A new model of wind erosion in the presence of vegetation
+    #Note that implementation only works in 1D currently
+
+    #Initialize shear variables and other grid parameters
+    ny, nx = xgrd.shape
+    ustar_okin = ustar.copy()
+
+    c1 = okin_params['okin_c1_veg']
+    intercept = okin_params['okin_initialred_veg']
+    if intercept > 1:
+        intercept = 1
+
+    #Calculate shear reduction by looking through all cells that have plants present and looking downwind of those features
+    for igridy in range(ny):
+        zp = hveg[igridy, :]
+        x = xgrd[igridy, :]
+        red_all = np.zeros(x.size)
+
+        for igrid in range(nx):
+
+            if zp[igrid] > 0:         # only look at cells with a roughness element
+                mult = np.ones(x.size)
+                red = np.zeros(x.size)
+                h = zp[igrid] #vegetation height at the appropriate cell
+
+                xrel = x - x[igrid]
+                #xrel = -xrel
+                for igrid2 in range(nx):
+
+                    if xrel[igrid2] >= 0:
+                        # apply okin model
+                        mult[igrid2] = intercept + (1 - intercept) * (1 - np.exp(-xrel[igrid2] * c1 / h))
+                    else:
+                        mult[igrid2] = 1
+
+                red = 1 - mult
+
+                # fix potential issues for summation
+                ix = red < 0.00001
+                red[ix] = 0
+                ix = red > 1
+                red[ix] = 1
+                ix = xrel < 0
+                red[ix] = 0
+
+                # combine all reductions between plants
+                red_all = red_all + red
+
+        # cant have more than 100% reduction
+        ifix = red_all > 1
+        red_all[ifix] = 1
+
+        #update shear velocity according to Okin (note does not operate on shear stress)
+        mult_all = 1 - red_all
+
+        ustarveg = ustar_okin[igridy, :] * mult_all
+
+        ix = ustarveg < 0.01
+        ustarveg[ix] = 0.01 #some small number so transport code doesnt crash
+
+        ustar_okin[igridy,:] = ustarveg
+
+    return ustar_okin
diff --git a/aeolis/wind.py b/aeolis/wind.py
index 45a6c1f6..0ea51ffc 100644
--- a/aeolis/wind.py
+++ b/aeolis/wind.py
@@ -33,6 +33,14 @@
 #import scipy.interpolate
 from scipy import ndimage, misc
 import matplotlib.pyplot as plt
+from scipy.integrate import quad
+from scipy.integrate import romberg
+from scipy.optimize import root_scalar
+from scipy.signal import convolve
+from scipy.special import sici
+from scipy.special import erfi
+from scipy.interpolate import griddata
+
 
 # package modules
 import aeolis.shear
@@ -73,10 +81,15 @@ def initialize(s, p):
     
     if p['process_shear']:
         if p['ny'] > 0:
-            s['shear'] = aeolis.shear.WindShear(s['x'], s['y'], s['zb'],
-                                                dx=p['dx'], dy=p['dy'],
-                                                L=p['L'], l=p['l'], z0=z0,
-                                                buffer_width=p['buffer_width'])
+            if p['method_shear'] == 'fft':
+                s['shear'] = aeolis.shear.WindShear(s['x'], s['y'], s['zb'],
+                                                    dx=p['dx'], dy=p['dy'],
+                                                    L=p['L'], l=p['l'], z0=z0,
+                                                    buffer_width=p['buffer_width'])
+            else:
+                s['shear'] = aeolis.rotation.rotationClass(s['x'], s['y'], s['zb'],
+                                                          dx=p['dx'], dy=p['dy'],
+                                                          buffer_width=100)
         else:
             s['shear'] = np.zeros(s['x'].shape)
 
@@ -84,7 +97,6 @@ def initialize(s, p):
 
 def interpolate(s, p, t):
     '''Interpolate wind velocity and direction to current time step
-
     Interpolates the wind time series for velocity and direction to
     the current time step. The cosine and sine of the direction angle
     are interpolated separately to prevent zero-crossing errors. The
@@ -220,8 +232,32 @@ def shear(s,p):
     
     # Compute shear velocity field (including separation)
 
-    if 'shear' in s.keys() and p['process_shear'] and p['ny'] > 0:
-        
+    if 'shear' in s.keys() and p['process_shear'] and p['ny'] > 0 and p['method_shear'] == 'duna2d':
+
+        shear_params = {'alfa': 3, 'beta': 1, 'c': p['c_b'], 'mu_b': p['mu_b'], 'sep_filter_iterations': p['sep_filter_iterations'], 'zsep_y_filter': p['zsep_y_filter'], 'process_separation': p['process_separation'], 'tau_sep': 0.5, 'slope': 0.2, 'rhoa': p['rhoa'], 'shear_type': p['method_shear']}
+
+        s['shear'](x=s['x'], y=s['y'], z=s['zb'], udir=s['udir'][0, 0], ustar=s['ustar'], tau=s['tau'],  hveg=s['hveg'], type_flag = 3, params=shear_params)
+
+        s['tau'] = s['shear'].get_tau()
+        s['taus'] = - s['tau'] * np.sin((-p['alfa'] + s['udir']) / 180. * np.pi)
+        s['taun'] = - s['tau'] * np.cos((-p['alfa'] + s['udir']) / 180. * np.pi)
+
+        s = stress_velocity(s, p)
+
+    elif 'shear' in s.keys() and p['process_shear'] and p['ny'] > 0 and p['method_shear'] == 'quasi2d':
+
+        z0 = calculate_z0(p, s)
+
+        shear_params = {'L': p['L'], 'l': p['l'], 'kappa': p['kappa'], 'viscNu': 1.5e-5, 'z0': z0, 'rho': p['rhoa'], 'sep_angle': 14, 'shear_type': p['method_shear'], 'process_separation': p['process_separation'], 'c': p['c_b'], 'mu_b': p['mu_b'], 'sep_filter_iterations': p['sep_filter_iterations'], 'zsep_y_filter': p['zsep_y_filter'], 'tau_sep': 0.5, 'slope': 0.2}
+
+        s['shear'](x=s['x'], y=s['y'], z=s['zb'], udir=s['udir'][0, 0], ustar=s['ustar'], tau=s['tau'],  hveg=s['hveg'], type_flag = 3, params=shear_params)
+        s['tau'] = s['shear'].get_tau()
+
+        s['taus'] = - s['tau'] * np.sin((-p['alfa'] + s['udir']) / 180. * np.pi)
+        s['taun'] = - s['tau'] * np.cos((-p['alfa'] + s['udir']) / 180. * np.pi)
+        s = stress_velocity(s, p)
+
+    elif 'shear' in s.keys() and p['process_shear'] and p['ny'] > 0:
         s['shear'](x=s['x'], y=s['y'], z=s['zb'],
                    taux=s['taus'], tauy=s['taun'],
                    u0=s['uw'][0,0], udir=s['udir'][0,0],
@@ -526,4 +562,459 @@ def separation1d(s, p):
     if direction == 2:
         zsepout = np.flip(zsepout)
 
-    return zsepout
\ No newline at end of file
+    return zsepout
+
+
+
+class Kroy2002():
+
+    def __init__(self,xgrid, z, z0, dune_length, rho, mu, tau, sep=False, sepInd=None, sep_ang=14, kappa=0.4):
+        self.xgrid = xgrid
+        self.z0 = z0
+        self.z = z
+        self.L = 0.5 * dune_length
+        self.Phi = self.L / z0
+        self.rho = rho
+        self.mu = mu
+        self._kappa = kappa
+        self.tau_inf = tau
+        self.tau_p = None
+        self.tau_tot = None
+        self.sep_slope = np.tan(np.radians(sep_ang))
+
+        self.dz = np.zeros(self.xgrid.shape)
+        ## check for equal spacing
+        dx = self.xgrid[1:] - self.xgrid[0:-1]
+        if np.all(np.abs(dx[0] - dx) <= 1e-6):
+            dx = dx[0]
+        else:
+            raise Exception('Unequally spaced dune profile data is not supported!')
+
+        ## use second order central differences except on ends
+        self.dz[0] = (-3.0 * self.z[0] + 4.0 * self.z[1] - self.z[2]) / (2.0 * dx)
+        self.dz[-1] = (3.0 * self.z[-1] - 4.0 * self.z[-2] + self.z[-3]) / (2.0 * dx)
+        self.dz[1:-1] = (self.z[2:] - self.z[0:-2]) / (2.0 * dx)
+
+        if sep and sepInd is None:
+            raise Exception('Must provide and index for the location of separation if "sep=True"')
+        self._sep = sep
+        self._sepInd = sepInd
+
+    def calc_tau_p(self):
+        # self._Phi = L/z0
+
+        # check for separation (max slope of -14 deg.)
+        #if self.dz.min() < -self.sep_slope or self._sep:
+        #    print('Separation possible based on dune profile!')
+        #    self._build_envelope()
+        #else:
+        #    self.env = self.z
+        self.env = self.z
+
+        # determine alpha and beta using eq. 13 from Kroy 03/2002
+        def phi_func(x):
+            # implicit equation for determining phi
+            val = x - (2.0 * self._kappa ** 2 * self.Phi / np.log(x))
+            deriv1 = 1 + (2.0 * self._kappa ** 2 * self.Phi / (x * np.log(x) ** 2))
+            # deriv2 = (2.0*self._kappa**2*Phi) * (np.log(x) + 2.0)/(x**2 * np.log(x)**3)
+            return val, deriv1  # , deriv2
+
+        phi_sol = root_scalar(phi_func,
+                              method='newton',
+                              x0=2.0 * self._kappa ** 2 * self.Phi / np.log(2.0 * self._kappa ** 2 * self.Phi),
+                              fprime=True)  # , fprime2=True)
+
+        if (phi_sol.converged):
+            phi = phi_sol.root
+        else:
+            exit('Implicit equation for phi did not converge! Exiting!')
+
+        fact = 1 + np.log(phi) + 2.0 * np.log(0.5 * np.pi) + 4.0 * np.euler_gamma
+        self._alpha = ((np.log(self.Phi ** 2 / np.log(self.Phi)) ** 2) / (2.0 * np.log(phi) ** 3)) * fact
+        self._beta = np.pi / fact
+
+        # calculate shear stress perturbation (2D) using eq. 16a from Kroy 03/2002
+        x = self.xgrid
+        L = self.L
+        fft_out = np.fft.rfft(self.env)
+        freq = 2.0 * np.pi * np.fft.rfftfreq(x.shape[0], x[1] - x[0])
+        fft_tau = self._alpha * (freq + 1j * self._beta * np.abs(freq)) * fft_out
+        self.tau_p = np.fft.irfft(fft_tau)
+        #self.tau_p_ana = 2.0 * self._alpha * (
+        #            np.pi ** (-0.5) - (x / L) * (self._beta + erfi(x / L)) * np.exp(-(x / L) ** 2))
+
+
+        #dx = x[1]-x[0]
+        #x2 = x[1:-1]-dx/2
+
+        #print(np.size(x))
+        #print(np.size(x2))
+        #print(np.size(self.tau_p))
+
+       # tau_p_interp = np.interp(x, x2, self.tau_p)
+
+        #func = np.interp1d(x2, self.tau_p,  axis=0,  # interpolate along columns
+        #                bounds_error=False,
+        #                kind='linear')
+
+        #tau_p_interp = func(x)
+
+        #self.tau_tot = self.tau_inf * (1 + tau_p_interp)
+
+    def calc_init_y(self, yPlus):
+        # calculate smallest needed y for desired y+
+        max_loc_u_star = np.sqrt(self.tau_tot.max() / self.rho)
+        first_y = (self.mu / self.rho) * yPlus / max_loc_u_star
+        return 2.0 * first_y  # y+ taken as centroid of cell, want to return total cell thickness
+
+    def S(self, x):
+        s = np.ones(x.shape)
+        s[x >= 0.5 * np.pi] = 0.0
+        s[x <= 0.5 * np.pi] = 0.0
+        return s
+
+    def _build_envelope(self):
+        def interp_end_vals(xr):
+            """Interpolates the elevation and slope at the reattachment point given by
+            xr"""
+
+            # get index of value just less than xr
+            #print(self.xgrid.max(), xr)
+            xrInd = np.argwhere(self.xgrid > xr)[0][0] - 1
+
+            # interp the values
+            fact = ((xr - self.xgrid[xrInd]) / (self.xgrid[xrInd + 1] - self.xgrid[xrInd]))
+            zr = self.z[xrInd] + fact * (self.z[xrInd + 1] - self.z[xrInd])
+            dzr = self.dz[xrInd] + fact * (self.dz[xrInd + 1] - self.dz[xrInd])
+
+            return zr, dzr
+
+        def solve_coef(xd, xr, zd, dzd, zr, dzr):
+            """Solves for the coefficients of a cubic polynomial given the value (zd, zr) and
+            derivative (dzd, dzr) at either end (xd, xr) """
+            # setup the system
+            A = np.array(
+                [[xd ** 3, xd ** 2, xd, 1.0],
+                 [3.0 * xd ** 2, 2.0 * xd, 1.0, 0.0],
+                 [xr ** 3, xr ** 2, xr, 1.0],
+                 [3.0 * xr ** 2, 2.0 * xr, 1.0, 0.0]]
+            )
+            b = [zd, dzd, zr, dzr]
+
+            # solve the system
+            return np.linalg.solve(A, b)
+
+        def sepLengthRoot(xr):
+            """Function for determining the separation bubble length which results in a
+            cubic polynomial with max slope equal to the limit slope"""
+            # interpolate endpoint values
+            zr, dzr = interp_end_vals(xr)
+
+            # solve for the polynomial
+            poly = solve_coef(xSepVal, xr, zd, dzd, zr, dzr)
+
+            # calc the max derivative
+            deriv_max = -(poly[1] ** 2 / (3.0 * poly[0])) + poly[2]
+
+            # return the difference between the deriv and the limit
+            return deriv_max - (-self.sep_slope)
+
+        # location of potential separation
+        if self._sepInd is None:
+            xSepInd = np.argwhere(self.dz < -self.sep_slope)[0][0] - 1
+        else:
+            xSepInd = self._sepInd
+
+        # useful values at separation point
+        xSepVal = self.xgrid[xSepInd]
+        zd = self.z[xSepInd]
+        # back-difference for left-side derivative at separation point
+        dzd = (self.z[xSepInd] - self.z[xSepInd - 1]) / (self.xgrid[xSepInd] - self.xgrid[xSepInd - 1])
+
+        # bubble length initial guess (Kroy2002 Eq. 29)
+        nu = dzd / self.sep_slope
+        Lb = 1.5 * (zd / self.sep_slope) * (
+                1.0 + 0.25 * nu + 0.125 * nu ** 2)
+
+        print('Predicted separation point: {:.4f}'.format(xSepVal))
+        print('Initial reattachment point: {:.4f}'.format(xSepVal + Lb))
+
+        # solve for actual reattachment length
+        # root_sol = root_scalar(sepLengthRoot,
+        #                        method='secant',
+        #                        x0=(xSepVal + Lb),
+        #                        x1=(xSepVal + Lb + 1e-6))
+        root_sol = root_scalar(sepLengthRoot,
+                               method='toms748',
+                               bracket=[xSepVal + 0.1, self.xgrid[-5]])
+        xr_fin = root_sol.root
+
+        print('\nSeparation Bubble Solve:')
+        print('\tInitial Separation Length: {:}'.format(Lb))
+        print('\tConverged: {:}'.format(root_sol.converged))
+        print('\tReason: {:}'.format(root_sol.flag))
+        print('\tSeparation Length: {:}'.format(xr_fin - xSepVal))
+        print('\tReattachment point: {:}'.format(xr_fin))
+
+        # fake it
+        # xr_fin = xSepVal + Lb
+        # solve for the final polynomial coefficients
+        zr, dzr = interp_end_vals(xr_fin)
+        poly = solve_coef(xSepVal, xr_fin, zd, dzd, zr, dzr)
+
+        # set the envelope
+        xrInd = np.argwhere(self.xgrid > xr_fin)[0][0] - 1
+        bub = lambda x: poly[0] * x ** 3 + poly[1] * x ** 2 + poly[2] * x + poly[3]
+
+        self.env = self.z.copy()
+        self.env[xSepInd:xrInd + 1] = bub(self.xgrid[xSepInd:xrInd + 1])
+
+def compute_shear_perturbation_1D(x, z, tau, params):
+    '''Compute wind shear perturbation for given free-flow wind
+    speed on computational grid. based on same implementation in Duna'''
+
+    ny, nx = x.shape
+    x = x[0,:]
+    dx = x[1] - x[0]
+    dx = np.abs(dx)
+    alfa = params['alfa']
+    beta = params['beta']
+    nx = nx-1
+    tau_update = tau.copy()
+
+
+    for igridy in range(ny):
+        zb = z[igridy, :]
+        dzbdx = np.zeros(x.size)
+
+        dzbdx[1:-1] = (zb[2:] - zb[0:-2]) / 2 / dx
+        tau_over_tau0 = np.zeros(x.size)
+
+        for i in range(nx + 1):
+            integ = 0
+            startval = i - nx
+            endval = i - 1
+            for j in np.arange(startval, endval + 1):
+                if j != 0:
+                    integ = integ + dzbdx[i - j] / (j * np.pi)
+            tau_over_tau0[i] = alfa * (integ + beta * dzbdx[i]) + 1
+            tau_over_tau0[i] = np.maximum(tau_over_tau0[i], 0.1)
+
+        tau_update[igridy,:] = tau[igridy,:] * tau_over_tau0
+
+    return tau_update
+
+def separation_quasi2d(x,y,z,tau,dx,dy, params):
+
+    if params['process_separation']:
+
+        c = params['c']
+        mu_b = params['mu_b']
+        sep_filter_iterations = params['sep_filter_iterations']
+        zsep_y_filter = params['zsep_y_filter']
+
+        # Initialize grid and bed dimensions
+        nx = len(z[1])
+        ny = len(z[0])
+        dx = dx
+        dy = dy
+        dzx = np.zeros(z.shape)
+        dzdx0 = np.zeros(z.shape)
+        dzdx1 = np.zeros(z.shape)
+        stall = np.zeros(z.shape)
+        bubble = np.zeros(z.shape)
+        k = np.array(range(0, nx))
+
+        zsep = np.zeros(z.shape)  # total separation bubble
+        zsep_new = np.zeros(z.shape)  # first-oder separation bubble surface
+
+        zfft = np.zeros((ny, nx), dtype=complex)
+
+        # Compute bed slope angle in x-dir
+        dzx[:, :-2] = np.rad2deg(np.arctan((z[:, 2:] - z[:, :-2]) / (2. * dx)))
+        dzx[:, -2] = dzx[:, -3]
+        dzx[:, -1] = dzx[:, -2]
+
+        # Determine location of separation bubbles
+        '''Separation bubble exist if bed slope angle (lee side) 
+        is larger than max angle that wind stream lines can 
+        follow behind an obstacle (mu_b = ..)'''
+
+        stall += np.logical_and(abs(dzx) > mu_b, dzx < 0.)
+
+        stall[:, 1:-1] += np.logical_and(stall[:, 1:-1] == 0, stall[:, :-2] > 0., stall[:, 2:] > 0.)
+
+        # Define separation bubble
+        bubble[:, :-1] = (stall[:, :-1] == 0.) * (stall[:, 1:] > 0.)
+
+        # Better solution for cleaner separation bubble, but no working Barchan dune (yet)
+        p = 1
+        bubble[:, p:] = bubble[:, :-p]
+        bubble[:, -p:] = 0
+
+        bubble = bubble.astype(int)
+
+        # Count separation bubbles
+        n = np.sum(bubble)
+        bubble_n = np.asarray(np.where(bubble == True)).T
+
+        # Walk through all separation bubbles and determine polynoms
+        j = 9999
+        for k in range(0, n):
+
+            i = bubble_n[k, 1]
+            j = bubble_n[k, 0]
+
+            # Bart: check for negative wind direction
+            #if np.sum(gc['taux']) >= 0:
+            idir = 1
+            #else:
+            #    idir = -1
+
+            ix_neg = (dzx[j, i + idir * 5:] >= 0)  # i + 5??
+
+            if np.sum(ix_neg) == 0:
+                zbrink = z[j, i]  # z level of brink at z(x0)
+            else:
+                zbrink = z[j, i] - z[j, i + idir * 5 + idir * np.where(ix_neg)[0][0]]
+
+            # Better solution and cleaner separation bubble, but no working Barchan dune (yet)
+            dzdx0 = (z[j, i] - z[j, i - 3]) / (3. * dx)
+
+            a = dzdx0 / c
+            ls = np.minimum(np.maximum((3. * zbrink / (2. * c) * (1. + a / 4. + a ** 2 / 8.)), 0.1), 200.)
+
+            a2 = -3 * zbrink / ls ** 2 - 2 * dzdx0 / ls
+            a3 = 2 * zbrink / ls ** 3 + dzdx0 / ls ** 2
+
+            i_max = min(i + int(ls / dx) + 1, int(nx - 1))
+
+            if idir == 1:
+                xs = x[j, i:i_max] - x[j, i]
+            else:
+                xs = -(x[j, i:i_max] - x[j, i])
+
+            zsep_new[j, i:i_max] = (a3 * xs ** 3 + a2 * xs ** 2 + dzdx0 * xs + z[j, i])
+
+            # Choose maximum of bedlevel, previous zseps and new zseps
+            zsep[j, :] = np.maximum.reduce([z[j, :], zsep[j, :], zsep_new[j, :]])
+
+            for filter_iter in range(sep_filter_iterations):
+
+                zsep_new = np.zeros(zsep.shape)
+
+                Cut = 1.5
+                dk = 2.0 * np.pi / (np.max(x))
+                zfft[j, :] = np.fft.fft(zsep[j, :])
+                zfft[j, :] *= np.exp(-(dk * k * dx) ** 2 / (2. * Cut ** 2))
+                zsep_fft = np.real(np.fft.ifft(zfft[j, :]))
+
+                if np.sum(ix_neg) == 0:
+                    zbrink = zsep_fft[i]
+                else:
+                    zbrink = zsep_fft[i] - zsep_fft[i + idir * 5 + idir * np.where(ix_neg)[0][0]]
+
+                # First order polynom
+                dzdx1 = (zsep_fft[i] - zsep_fft[i - 3]) / (3. * dx)
+
+                a = dzdx1 / c
+
+                ls = np.minimum(np.maximum((3. * zbrink / (2. * c) * (1. + a / 4. + a ** 2 / 8.)), 0.1), 200.)
+
+                a2 = -3 * zbrink / ls ** 2 - 2 * dzdx1 / ls
+                a3 = 2 * zbrink / ls ** 3 + dzdx1 / ls ** 2
+
+                i_max1 = min(i + idir * int(ls / dx), int(nx - 1))
+
+                if idir == 1:
+                    xs1 = x[j, i:i_max1] - x[j, i]
+                else:
+                    xs1 = -(x[j, i:i_max1] - x[j, i])
+
+                zsep_new[j, i:i_max1] = (a3 * xs1 ** 3 + a2 * xs1 ** 2 + dzdx1 * xs1 + zbrink)
+
+                # Pick the maximum seperation bubble hieght at all locations
+                zsep[j, :] = np.maximum.reduce([z[j, :], zsep[j, :], zsep_new[j, :]])
+
+        # Smooth surface of separation bubbles over y direction
+        if zsep_y_filter:
+            zsep = ndimage.gaussian_filter1d(zsep, sigma=0.2, axis=0)
+
+        # Correct for any seperation bubbles that are below the bed surface following smoothing
+        ilow = zsep < z
+        zsep[ilow] = z[ilow]
+
+        hsep = zsep - z
+        tau_sep = params['tau_sep']
+        slope = params['slope']
+        delta = 1. / (slope * tau_sep)
+
+        zsepdelta = np.minimum(np.maximum(1. - delta * hsep, 0.), 1.)
+
+        tau = tau * zsepdelta
+
+    else:
+
+        zsep = np.zeros(z.shape)
+        hsep = np.zeros(z.shape)
+        tau = tau
+
+    return zsep, hsep, tau
+
+def compute_shear_perturbation_kroy1D(x, y, z, tau, params):
+    '''Compute wind shear perturbation for given free-flow wind
+    speed on computational grid. based on same implementation in Duna'''
+
+    ny, nx = x.shape
+    x = x[0,:]
+    dx = x[1] - x[0]
+    dx = np.abs(dx)
+
+    x2 = x[1:]-dx/2
+
+    L = params['L']
+    l = params['l']
+    kappa = params['kappa']
+    mu = params['rho']*params['viscNu']
+    rho = params['rho']
+    sep_angle = params['sep_angle'] #14
+    z0 = params['z0']
+
+    tau_update = tau.copy()
+    ix = np.isnan(tau_update)
+    tau_update[ix] = 0
+
+    for igridy in range(ny):
+        zb = z[igridy, :]
+        tau_in = tau[igridy, :]
+        kroy = Kroy2002(x, zb, z0, L, rho, mu, tau_in, sep=False, sepInd=None, sep_ang=sep_angle, kappa=kappa)
+        kroy.calc_tau_p()
+        #tau_kroy = kroy.tau_tot
+        tau_p_kroy = kroy.tau_p
+        try:
+            tau_p_interp = np.interp(x, x2, tau_p_kroy)
+        except:
+            tau_p_interp = tau_p_kroy
+
+        ix = np.isnan(tau_p_interp)
+        tau_p_interp[ix] = 0
+        tau_kroy = tau_in * (1 + tau_p_interp)
+
+        ix = np.isinf(tau_kroy)
+        tau_kroy[ix] = tau_in[ix]
+
+        tau_update[igridy, :] = tau_kroy
+
+    #fix bad values
+    ix = np.isinf(tau_update)
+    tau_update[ix] = tau[ix]
+    ix = np.isnan(tau_update)
+    tau_update[ix] = tau[ix]
+    ix = tau_update < 0
+    tau_update[ix] = tau[ix]
+
+    xdata = x
+    ydata = y
+    tau_data = tau_update
+
+    return tau_update
\ No newline at end of file