diff --git a/wetb/hawc2/htc_contents.py b/wetb/hawc2/htc_contents.py
index 4e184974dfdd626c0099d8ff065f6f68d2bf95f1..a38b08ec515895d3fc2bfd2de91b073410e9f76b 100644
--- a/wetb/hawc2/htc_contents.py
+++ b/wetb/hawc2/htc_contents.py
@@ -326,130 +326,4 @@ class HTCSensor(HTCLine):
                                 ("", "\t" + self.str_values())[bool(self.values)],
                                 ("", "\t" + self.comments)[bool(self.comments.strip())])
 
-class HTCDefaults(object):
-
-
-    empty_htc = """begin simulation;
-        time_stop 600;
-        solvertype    1;    (newmark)
-        on_no_convergence continue;
-        convergence_limits 1E3 1.0 1E-7; ; . to run again, changed 07/11
-        begin newmark;
-          deltat    0.02;
-        end newmark;
-    end simulation;
-    ;
-    ;----------------------------------------------------------------------------------------------------------------------------------------------------------------
-    ;
-    begin new_htc_structure;
-      begin orientation;
-      end orientation;
-      begin constraint;
-      end constraint;
-    end new_htc_structure;
-    ;
-    ;----------------------------------------------------------------------------------------------------------------------------------------------------------------
-    ;
-    begin wind ;
-      density                 1.225 ;
-      wsp                     10   ;
-      tint                    1;
-      horizontal_input        1     ;            0=false, 1=true
-      windfield_rotations     0 0.0 0.0 ;    yaw, tilt, rotation
-      center_pos0             0 0 -30 ; hub heigth
-      shear_format            1   0;0=none,1=constant,2=log,3=power,4=linear
-      turb_format             0     ;  0=none, 1=mann,2=flex
-      tower_shadow_method     0     ;  0=none, 1=potential flow, 2=jet
-    end wind;
-    ;
-    ;----------------------------------------------------------------------------------------------------------------------------------------------------------------
-    ;
-    begin dll;
-    end dll;
-    ;
-    ;----------------------------------------------------------------------------------------------------------------------------------------------------------------
-    ;
-    begin output;
-      general time;
-    end output;
-    exit;"""
-
-
-    def add_mann_turbulence(self, L=29.4, ae23=1, Gamma=3.9, seed=1001, high_frq_compensation=True,
-                            filenames=None,
-                            no_grid_points=(16384, 32, 32), box_dimension=(6000, 100, 100),
-                            dont_scale=False,
-                            std_scaling=None):
-        wind = self.add_section('wind')
-        wind.turb_format = 1
-        mann = wind.add_section('mann')
-        if 'create_turb_parameters' in mann:
-            mann.create_turb_parameters.values = [L, ae23, Gamma, seed, int(high_frq_compensation)]
-        else:
-            mann.add_line('create_turb_parameters', [L, ae23, Gamma, seed, int(high_frq_compensation)], "L, alfaeps, gamma, seed, highfrq compensation")
-        if filenames is None:
-            fmt = "mann_l%.1f_ae%.2f_g%.1f_h%d_%dx%dx%d_%.3fx%.2fx%.2f_s%04d%c.turb"
-            import numpy as np
-            dxyz = tuple(np.array(box_dimension) / no_grid_points)
-            filenames = ["./turb/" + fmt % ((L, ae23, Gamma, high_frq_compensation) + no_grid_points + dxyz + (seed, uvw)) for uvw in ['u', 'v', 'w']]
-        if isinstance(filenames, str):
-            filenames = ["./turb/%s_s%04d%s.bin" % (filenames, seed, c) for c in ['u', 'v', 'w']]
-        for filename, c in zip(filenames, ['u', 'v', 'w']):
-            setattr(mann, 'filename_%s' % c, filename)
-        for c, n, dim in zip(['u', 'v', 'w'], no_grid_points, box_dimension):
-            setattr(mann, 'box_dim_%s' % c, "%d %.4f" % (n, dim / (n - 1)))
-        if dont_scale:
-            mann.dont_scale = 1
-        else:
-            try:
-                del mann.dont_scale
-            except KeyError:
-                pass
-        if std_scaling is not None:
-            mann.std_scaling = "%f %f %f" % std_scaling
-        else:
-            try:
-                del mann.std_scaling
-            except KeyError:
-                pass
-            
-
-
-    def add_turb_export(self, filename="export_%s.turb", samplefrq = None):
-        exp = self.wind.add_section('turb_export', allow_duplicate=True)
-        for uvw in 'uvw':
-            exp.add_line('filename_%s'%uvw, [filename%uvw])
-        sf = samplefrq or max(1,int( self.wind.mann.box_dim_u[1]/(self.wind.wsp[0] * self.deltat())))
-        exp.samplefrq = sf
-        if "time" in self.output:
-            exp.time_start = self.output.time[0]
-        else:
-            exp.time_start = 0
-        exp.nsteps = (self.simulation.time_stop[0]-exp.time_start[0]) / self.deltat()
-        for vw in 'vw':
-            exp.add_line('box_dim_%s'%vw, self.wind.mann['box_dim_%s'%vw].values)
-
-        
-
-    def import_dtu_we_controller_input(self, filename):
-        dtu_we_controller = [dll for dll in self.dll if dll.name[0] == 'dtu_we_controller'][0]
-        with open (filename) as fid:
-            lines = fid.readlines()
-        K_r1 = float(lines[1].replace("K = ", '').replace("[Nm/(rad/s)^2]", ''))
-        Kp_r2 = float(lines[4].replace("Kp = ", '').replace("[Nm/(rad/s)]", ''))
-        Ki_r2 = float(lines[5].replace("Ki = ", '').replace("[Nm/rad]", ''))
-        Kp_r3 = float(lines[7].replace("Kp = ", '').replace("[rad/(rad/s)]", ''))
-        Ki_r3 = float(lines[8].replace("Ki = ", '').replace("[rad/rad]", ''))
-        KK = lines[9].split("]")
-        KK1 = float(KK[0].replace("K1 = ", '').replace("[deg", ''))
-        KK2 = float(KK[1].replace(", K2 = ", '').replace("[deg^2", ''))
-        cs = dtu_we_controller.init
-        cs.constant__11.values[1] = "%.6E" % K_r1
-        cs.constant__12.values[1] = "%.6E" % Kp_r2
-        cs.constant__13.values[1] = "%.6E" % Ki_r2
-        cs.constant__16.values[1] = "%.6E" % Kp_r3
-        cs.constant__17.values[1] = "%.6E" % Ki_r3
-        cs.constant__21.values[1] = "%.6E" % KK1
-        cs.constant__22.values[1] = "%.6E" % KK2
-
 
diff --git a/wetb/hawc2/htc_extensions.py b/wetb/hawc2/htc_extensions.py
new file mode 100644
index 0000000000000000000000000000000000000000..67b9f7fe0bccbc5b08d9dfb1268cc9eb2058355e
--- /dev/null
+++ b/wetb/hawc2/htc_extensions.py
@@ -0,0 +1,160 @@
+'''
+Created on 20/01/2014
+
+@author: MMPE
+
+See documentation of HTCFile below
+
+'''
+from __future__ import division
+from __future__ import unicode_literals
+from __future__ import print_function
+from __future__ import absolute_import
+from builtins import zip
+from builtins import int
+from builtins import str
+from future import standard_library
+import os
+from wetb.wind.shear import log_shear, power_shear
+standard_library.install_aliases()
+
+
+
+class HTCDefaults(object):
+
+
+    empty_htc = """begin simulation;
+        time_stop 600;
+        solvertype    1;    (newmark)
+        on_no_convergence continue;
+        convergence_limits 1E3 1.0 1E-7; ; . to run again, changed 07/11
+        begin newmark;
+          deltat    0.02;
+        end newmark;
+    end simulation;
+    ;
+    ;----------------------------------------------------------------------------------------------------------------------------------------------------------------
+    ;
+    begin new_htc_structure;
+      begin orientation;
+      end orientation;
+      begin constraint;
+      end constraint;
+    end new_htc_structure;
+    ;
+    ;----------------------------------------------------------------------------------------------------------------------------------------------------------------
+    ;
+    begin wind ;
+      density                 1.225 ;
+      wsp                     10   ;
+      tint                    1;
+      horizontal_input        1     ;            0=false, 1=true
+      windfield_rotations     0 0.0 0.0 ;    yaw, tilt, rotation
+      center_pos0             0 0 -30 ; hub heigth
+      shear_format            1   0;0=none,1=constant,2=log,3=power,4=linear
+      turb_format             0     ;  0=none, 1=mann,2=flex
+      tower_shadow_method     0     ;  0=none, 1=potential flow, 2=jet
+    end wind;
+    ;
+    ;----------------------------------------------------------------------------------------------------------------------------------------------------------------
+    ;
+    begin dll;
+    end dll;
+    ;
+    ;----------------------------------------------------------------------------------------------------------------------------------------------------------------
+    ;
+    begin output;
+      general time;
+    end output;
+    exit;"""
+
+
+    def add_mann_turbulence(self, L=29.4, ae23=1, Gamma=3.9, seed=1001, high_frq_compensation=True,
+                            filenames=None,
+                            no_grid_points=(16384, 32, 32), box_dimension=(6000, 100, 100),
+                            dont_scale=False,
+                            std_scaling=None):
+        wind = self.add_section('wind')
+        wind.turb_format = 1
+        mann = wind.add_section('mann')
+        if 'create_turb_parameters' in mann:
+            mann.create_turb_parameters.values = [L, ae23, Gamma, seed, int(high_frq_compensation)]
+        else:
+            mann.add_line('create_turb_parameters', [L, ae23, Gamma, seed, int(high_frq_compensation)], "L, alfaeps, gamma, seed, highfrq compensation")
+        if filenames is None:
+            fmt = "mann_l%.1f_ae%.2f_g%.1f_h%d_%dx%dx%d_%.3fx%.2fx%.2f_s%04d%c.turb"
+            import numpy as np
+            dxyz = tuple(np.array(box_dimension) / no_grid_points)
+            filenames = ["./turb/" + fmt % ((L, ae23, Gamma, high_frq_compensation) + no_grid_points + dxyz + (seed, uvw)) for uvw in ['u', 'v', 'w']]
+        if isinstance(filenames, str):
+            filenames = ["./turb/%s_s%04d%s.bin" % (filenames, seed, c) for c in ['u', 'v', 'w']]
+        for filename, c in zip(filenames, ['u', 'v', 'w']):
+            setattr(mann, 'filename_%s' % c, filename)
+        for c, n, dim in zip(['u', 'v', 'w'], no_grid_points, box_dimension):
+            setattr(mann, 'box_dim_%s' % c, "%d %.4f" % (n, dim / (n - 1)))
+        if dont_scale:
+            mann.dont_scale = 1
+        else:
+            try:
+                del mann.dont_scale
+            except KeyError:
+                pass
+        if std_scaling is not None:
+            mann.std_scaling = "%f %f %f" % std_scaling
+        else:
+            try:
+                del mann.std_scaling
+            except KeyError:
+                pass
+            
+
+
+    def add_turb_export(self, filename="export_%s.turb", samplefrq = None):
+        exp = self.wind.add_section('turb_export', allow_duplicate=True)
+        for uvw in 'uvw':
+            exp.add_line('filename_%s'%uvw, [filename%uvw])
+        sf = samplefrq or max(1,int( self.wind.mann.box_dim_u[1]/(self.wind.wsp[0] * self.deltat())))
+        exp.samplefrq = sf
+        if "time" in self.output:
+            exp.time_start = self.output.time[0]
+        else:
+            exp.time_start = 0
+        exp.nsteps = (self.simulation.time_stop[0]-exp.time_start[0]) / self.deltat()
+        for vw in 'vw':
+            exp.add_line('box_dim_%s'%vw, self.wind.mann['box_dim_%s'%vw].values)
+
+        
+
+    def import_dtu_we_controller_input(self, filename):
+        dtu_we_controller = [dll for dll in self.dll if dll.name[0] == 'dtu_we_controller'][0]
+        with open (filename) as fid:
+            lines = fid.readlines()
+        K_r1 = float(lines[1].replace("K = ", '').replace("[Nm/(rad/s)^2]", ''))
+        Kp_r2 = float(lines[4].replace("Kp = ", '').replace("[Nm/(rad/s)]", ''))
+        Ki_r2 = float(lines[5].replace("Ki = ", '').replace("[Nm/rad]", ''))
+        Kp_r3 = float(lines[7].replace("Kp = ", '').replace("[rad/(rad/s)]", ''))
+        Ki_r3 = float(lines[8].replace("Ki = ", '').replace("[rad/rad]", ''))
+        KK = lines[9].split("]")
+        KK1 = float(KK[0].replace("K1 = ", '').replace("[deg", ''))
+        KK2 = float(KK[1].replace(", K2 = ", '').replace("[deg^2", ''))
+        cs = dtu_we_controller.init
+        cs.constant__11.values[1] = "%.6E" % K_r1
+        cs.constant__12.values[1] = "%.6E" % Kp_r2
+        cs.constant__13.values[1] = "%.6E" % Ki_r2
+        cs.constant__16.values[1] = "%.6E" % Kp_r3
+        cs.constant__17.values[1] = "%.6E" % Ki_r3
+        cs.constant__21.values[1] = "%.6E" % KK1
+        cs.constant__22.values[1] = "%.6E" % KK2
+
+
+class HTCExtensions(object):
+    def get_shear(self):
+        shear_type, parameter = self.wind.shear_format.values
+        z0 = -self.wind.center_pos0[2]
+        wsp = self.wind.wsp[0]
+        if shear_type==1: #constant
+            return lambda z : parameter
+        elif shear_type==3:
+            return power_shear(parameter, z0, wsp)
+        else:
+            raise NotImplementedError
\ No newline at end of file
diff --git a/wetb/hawc2/htc_file.py b/wetb/hawc2/htc_file.py
index 7191975254a2d4aa5bcf5db6fa6cad437caf7ba7..f783408e410388868e65a9dff5b4405f27139872 100644
--- a/wetb/hawc2/htc_file.py
+++ b/wetb/hawc2/htc_file.py
@@ -18,8 +18,8 @@ from wetb.utils.cluster_tools.cluster_resource import unix_path_old
 standard_library.install_aliases()
 from collections import OrderedDict
 
-from wetb.hawc2.htc_contents import HTCContents, HTCSection, HTCLine, \
-    HTCDefaults
+from wetb.hawc2.htc_contents import HTCContents, HTCSection, HTCLine
+from wetb.hawc2.htc_extensions import HTCDefaults, HTCExtensions
 import os
 from copy import copy
 
@@ -27,7 +27,7 @@ from copy import copy
 def fmt_path(path):
     return path.lower().replace("\\","/")
 
-class HTCFile(HTCContents, HTCDefaults):
+class HTCFile(HTCContents, HTCDefaults, HTCExtensions):
     """Wrapper for HTC files
 
     Examples:
diff --git a/wetb/hawc2/shear_file.py b/wetb/hawc2/shear_file.py
index 7c82cb14dac56008aa6db655581b5caa3acf0802..11d6d60803f5eb35f6d83bc231381cee8e2e91eb 100644
--- a/wetb/hawc2/shear_file.py
+++ b/wetb/hawc2/shear_file.py
@@ -10,62 +10,168 @@ from __future__ import absolute_import
 from builtins import range
 from io import open
 from future import standard_library
+from wetb.hawc2.htc_file import HTCFile
 standard_library.install_aliases()
-
 import numpy as np
 import os
 
-def save(filename, x_coordinates, z_coordinates, u=None, v=None, w=None):
-    """
-    Parameters
-    ----------
-    filename : str
-        filename
-    x_coordinates : array_like
-        lateral coordinates
-    z_coordinates : array_like
-        vertical coordinates
-    u : array_like, optional
-        shear_u component, normalized with U_mean\n
-        shape must be (#z_coordinates, #x_coordinates) or (#z_coordinates,)
-    v : array_like, optional
-        shear_v component, normalized with U_mean\n
-        shape must be (#z_coordinates, #x_coordinates) or (#z_coordinates,)
-    w : array_like, optional
-        shear_w component, normalized with U_mean\n
-        shape must be (#z_coordinates, #x_coordinates) or (#z_coordinates,)
-    """
 
-    shape = (len(z_coordinates), len(x_coordinates))
-    vuw = [v, u, w]
-    for i in range(3):
-        if vuw[i] is None:
-            if i == 1:
-                vuw[i] = np.ones((shape))
+class ShearFile(object):
+    """HAWC2 user defined shear file
+    
+    Examples:
+    ---------
+    
+    >>> sf = ShearFile([-55, 55], [30, 100, 160] , u=np.array([[0.7, 1, 1.3], [0.7, 1, 1.3]]).T)
+    >>> print (sf.uvw([-55,0],[65,135])) #uvw factors
+    [array([ 0.85 ,  1.175]), array([ 0.,  0.]), array([ 0.,  0.])]
+    >>> from wetb.wind.shear import power_shear
+    >>> print (sf.uvw([-55,0],[65,135], shear=power_shear(.2,100,10))) # uvw wind speeds
+    [array([  7.79832978,  12.47684042]), array([ 0.,  0.]), array([ 0.,  0.])]
+    >>> sf.save('test.dat')
+    """
+    
+    def __init__(self,v_positions, w_positions, u=None, v=None, w=None, shear=None):
+        """
+        Parameters
+        ----------
+        v_positions : array_like
+            lateral coordinates
+        w_positions : array_like
+            vertical coordinates
+        u : array_like, optional
+            shear_u component, normalized with U_mean\n
+            shape must be (#w_positions, #v_positions) or (#w_positions,)
+        v : array_like, optional
+            shear_v component, normalized with U_mean\n
+            shape must be (#w_positions, #v_positions) or (#w_positions,)
+        w : array_like, optional
+            shear_w component, normalized with U_mean\n
+            shape must be (#w_positions, #v_positions) or (#w_positions,)
+        """
+        self.v_positions = v_positions
+        self.w_positions = w_positions
+        shape = (len(w_positions), len(v_positions))
+        uvw = [u, v, w]
+        for i in range(3):
+            if uvw[i] is None:
+                if i == 0:
+                    uvw[i] = np.ones((shape))
+                else:
+                    uvw[i] = np.zeros((shape))
             else:
-                vuw[i] = np.zeros((shape))
-        else:
-            vuw[i] = np.array(vuw[i])
-            if len(vuw[i].shape) == 1 and vuw[i].shape[0] == shape[0]:
-                vuw[i] = np.repeat(np.atleast_2d(vuw[i]).T, shape[1], 1)
+                uvw[i] = np.array(uvw[i])
+                if len(uvw[i].shape) == 1 and uvw[i].shape[0] == shape[0]:
+                    uvw[i] = np.repeat(np.atleast_2d(uvw[i]).T, shape[1], 1)
+    
+                assert uvw[i].shape == shape, (i, uvw[i].shape, shape)
+        self.u, self.v, self.w = uvw
+        self.shear = shear
+    
+    def uvw(self, v,w,shear=None):
+        """Calculate u,v,w wind speeds at position(s) (v,w)
+        Parameters
+        ----------
+        v : int, float or array_like
+            v-coordinate(s)
+        w : int, float or array_like
+            w-coordinates(s)
+        shear : function or None
+            if function: f(height)->wsp
+            if None: self.shear is used, if not None. 
+            Otherwise wind speed factors instead of absolute wind speeds are returned
+            
+        Returns
+        -------
+        u,v,w 
+            wind speed(s) or wind speed factor(s) if shear not defined
+        """
+        shear = shear or self.shear or (lambda z: 1)
+        from scipy.interpolate import RegularGridInterpolator
+        wv = np.array([w,v]).T
+        return [RegularGridInterpolator((self.w_positions, self.v_positions), uvw)(wv)*shear(w) 
+                            for uvw in [self.u, self.v, self.w] if uvw is not None]
+        
+
+
+    def save(self, filename):
+        """Save user defined shear file
+        Parameters
+        ----------
+        filename : str:
+            Filename
+        """
+        
+        # exist_ok does not exist in Python27
+        filename = os.path.abspath(filename)
+        if not os.path.exists(os.path.dirname(filename)):
+            os.makedirs(os.path.dirname(filename))#, exist_ok=True)
+        with open(filename, 'w', encoding='utf-8') as fid:
+            fid.write(" # autogenerated shear file\n")
+            fid.write("  %d %d\n" % (len(self.v_positions), len(self.w_positions)))
+            
+            for i, (l,vuw) in enumerate(zip(['v', 'u', 'w'],[self.v,self.u,self.w])):
+                fid.write(" # shear %s component\n  " % l)
+                fid.write("\n  ".join([" ".join(["%.10f" % v for v in r ]) for r in vuw]))
+                fid.write("\n")
+            for yz, coor in (['v', self.v_positions], ['w', self.w_positions]):
+                fid.write(" # %s coordinates\n  " % yz)
+                fid.write("\n  ".join("%.10f" % v for v in coor))
+                fid.write("\n")
+
+    @staticmethod
+    def load(filename):
+        """Load shear file
+        Parameters 
+        ----------
+        filename : str
+            Filename
+            
+        Returns
+        -------
+        shear file : ShearFile-object
+        """
+        with open(filename) as fid:
+            lines = fid.readlines()
+        no_V,no_W=map(int,lines[1].split())
+        v,u,w = [np.array([row.split() for row in lines[3+(no_W+1)*i:3+(no_W+1)*(i+1)-1] ],dtype=float) for i in range(3)]
+        v_positions = np.array(lines[3+(no_W+1)*3:3+(no_W+1)*3+no_V],dtype=float)
+        w_positions = np.array(lines[3+(no_W+1)*3+(no_V+1):3+(no_W+1)*3+(no_V+1)+no_W],dtype=float)
+        return ShearFile(v_positions,w_positions,u,v,w)
+
+    @staticmethod
+    def load_from_htc(htc_file):
+        """Load shear file from HTC file including shear function
+        Parameters 
+        ----------
+        htc_file : str or HTCFile
+            Filename or HTCFile
+            
+        Returns
+        -------
+        shear file : ShearFile-object
+        """
+        if isinstance(htc_file,str):
+            htc_file = HTCFile(htc_file)
+        user_defined_shear_filename = os.path.join(htc_file.modelpath, htc_file.wind.user_defined_shear[0])
+        shear_file = ShearFile.load(user_defined_shear_filename)
+        shear_file.shear = htc_file.get_shear()
+        return shear_file
+
+def save(filename, v_coordinates, w_coordinates, u=None, v=None, w=None):
+    """Save shear file (deprecated)"""
+    ShearFile(v_coordinates, w_coordinates,u,v,w).save(filename)
+
 
-            assert vuw[i].shape == shape, (i, vuw[i].shape, shape)
 
-    # exist_ok does not exist in Python27
-    if not os.path.exists(os.path.dirname(filename)):
-        os.makedirs(os.path.dirname(filename))#, exist_ok=True)
-    with open(filename, 'w', encoding='utf-8') as fid:
-        fid.write(" # autogenerated shear file\n")
-        fid.write("  %d %d\n" % (len(x_coordinates), len(z_coordinates)))
-        for i, l in enumerate(['v', 'u', 'w']):
-            fid.write(" # shear %s component\n  " % l)
-            fid.write("\n  ".join([" ".join(["%.10f" % v for v in r ]) for r in vuw[i]]))
-            fid.write("\n")
-        for yz, coor in (['v', x_coordinates], ['w', z_coordinates]):
-            fid.write(" # %s coordinates\n  " % yz)
-            fid.write("\n  ".join("%.10f" % v for v in coor))
-            fid.write("\n")
 
 
 if __name__ == "__main__":
-    save("test.dat", [-55, 55], [30, 100, 160] , u=np.array([[0.7, 1, 1.3], [0.7, 1, 1.3]]).T)
+    from wetb.wind.shear import power_shear
+    sf = ShearFile([-55, 55], [30, 100, 160] , u=np.array([[0.7, 1, 1.3], [0.7, 1, 1.3]]).T)
+    print (sf.uvw([-55,0],[65,135])) #uvw factors
+    #[array([ 0.85 ,  1.175]), array([ 0.,  0.]), array([ 0.,  0.])]
+    print (sf.uvw([-55,0],[65,135], shear=power_shear(.2,100,10))) # uvw wind speeds
+    #[array([  7.79832978,  12.47684042]), array([ 0.,  0.]), array([ 0.,  0.])]
+    sf.save('test.dat')
+    
diff --git a/wetb/hawc2/tests/test_files/data/user_shear.dat b/wetb/hawc2/tests/test_files/data/user_shear.dat
new file mode 100644
index 0000000000000000000000000000000000000000..d8371b679430d84a5273bc6356871f24cb218c44
--- /dev/null
+++ b/wetb/hawc2/tests/test_files/data/user_shear.dat
@@ -0,0 +1,21 @@
+ # autogenerated shear file
+  2 3
+ # shear v component
+  0.0000000000 0.0000000000
+  0.0000000000 0.0000000000
+  0.0000000000 0.0000000000
+ # shear u component
+  0.8000000000 0.6000000000
+  1.0000000000 1.0000000000
+  1.2000000000 1.4000000000
+ # shear w component
+  0.0000000000 0.0000000000
+  0.0000000000 0.0000000000
+  0.0000000000 0.0000000000
+ # v coordinates
+  -55.0000000000
+  55.0000000000
+ # w coordinates
+  30.0000000000
+  100.0000000000
+  160.0000000000
diff --git a/wetb/hawc2/tests/test_files/htcfiles/test.htc b/wetb/hawc2/tests/test_files/htcfiles/test.htc
index 2dfe96ecf93e5d8a583461ed76a2c66524f37b96..ef87ececc5d83d5f98a8e43295216e685ba9673b 100644
--- a/wetb/hawc2/tests/test_files/htcfiles/test.htc
+++ b/wetb/hawc2/tests/test_files/htcfiles/test.htc
@@ -260,6 +260,7 @@ begin wind;
   tower_shadow_method     3;  0=none, 1=potential flow, 2=jet
   scale_time_start        0; 
   wind_ramp_factor   0.0 100 0.8 1.0;
+  user_defined_shear      ./data/user_shear.dat;
 ; iec_gust;
 ;
 ;  wind_ramp_abs  400.0  401.0  0.0   1.0;   wsp. after the step:  5.0 
diff --git a/wetb/hawc2/tests/test_htc_extensions.py b/wetb/hawc2/tests/test_htc_extensions.py
new file mode 100644
index 0000000000000000000000000000000000000000..fe84fbc2a1f4fb44dfeede46c4ed3adbe5c4c82b
--- /dev/null
+++ b/wetb/hawc2/tests/test_htc_extensions.py
@@ -0,0 +1,33 @@
+'''
+Created on 17/07/2014
+
+@author: MMPE
+'''
+from __future__ import print_function
+from __future__ import unicode_literals
+from __future__ import division
+from __future__ import absolute_import
+from io import open
+from builtins import str
+from builtins import zip
+from future import standard_library
+standard_library.install_aliases()
+import os
+import unittest
+
+from datetime import datetime
+from wetb.hawc2.htc_file import HTCFile, HTCLine
+
+import numpy as np
+
+
+tfp = os.path.join(os.path.dirname(__file__), 'test_files/htcfiles/')  # test file path
+class TestHtcFile(unittest.TestCase):
+
+    def test_get_shear(self):
+        htc = HTCFile(tfp+'test.htc')
+        self.assertEqual(htc.get_shear()(100), 10*(100/119)**.2)
+    
+if __name__ == "__main__":
+    #import sys;sys.argv = ['', 'Test.testName']
+    unittest.main()
diff --git a/wetb/hawc2/tests/test_htc_file.py b/wetb/hawc2/tests/test_htc_file.py
index 6b9caa9c7774ee2d9fb50c7f9776b41c638541c7..79b089272d8b8b2f73a8f8f2428c64101b024646 100644
--- a/wetb/hawc2/tests/test_htc_file.py
+++ b/wetb/hawc2/tests/test_htc_file.py
@@ -244,6 +244,7 @@ end turb_export;"""
                   './control/mech_brake.dll',
                   './control/servo_with_limits.dll',
                   './control/towclearsens.dll',
+                  './data/user_shear.dat',
                   self.testfilepath.replace("\\","/") + 'test.htc'
                   ]:
             try:
@@ -313,7 +314,6 @@ end turb_export;"""
         htc = HTCFile(self.testfilepath + "test_2xoutput.htc","../")
         self.assertEqual(len(htc.res_file_lst()), 4)
         
-    
 
 if __name__ == "__main__":
     #import sys;sys.argv = ['', 'Test.testName']
diff --git a/wetb/hawc2/tests/test_shear_file.py b/wetb/hawc2/tests/test_shear_file.py
index ade7c19a74c50d95dd44431d589253f07795f553..302cae18134dac29b34ac36c4232b6cbb02de9eb 100644
--- a/wetb/hawc2/tests/test_shear_file.py
+++ b/wetb/hawc2/tests/test_shear_file.py
@@ -9,18 +9,19 @@ from __future__ import division
 from __future__ import absolute_import
 from io import open
 from future import standard_library
+from wetb.hawc2.shear_file import ShearFile
 standard_library.install_aliases()
 import unittest
 from wetb.hawc2 import shear_file
 import numpy as np
 import os
 import shutil
-testfilepath = 'test_files/'
+tfp = os.path.join(os.path.dirname(__file__), 'test_files/')
 class TestShearFile(unittest.TestCase):
 
 
-    def test_shearfile(self):
-        f = testfilepath + "tmp_shearfile1.dat"
+    def test_shearfile_save(self):
+        f = tfp + "tmp_shearfile1.dat"
         shear_file.save(f, [-55, 55], [30, 100, 160] , u=np.array([[0.7, 1, 1.3], [0.7, 1, 1.3]]).T)
         with open(f) as fid:
             self.assertEqual(fid.read(),
@@ -50,7 +51,7 @@ class TestShearFile(unittest.TestCase):
 
 
     def test_shearfile2(self):
-        f = testfilepath + "tmp_shearfile2.dat"
+        f = tfp + "tmp_shearfile2.dat"
         shear_file.save(f, [-55, 55], [30, 100, 160] , u=np.array([0.7, 1, 1.3]).T)
         with open(f) as fid:
             self.assertEqual(fid.read(),
@@ -79,10 +80,20 @@ class TestShearFile(unittest.TestCase):
         os.remove(f)
 
     def test_shear_makedirs(self):
-        f = testfilepath + "shear/tmp_shearfile2.dat"
+        f = tfp + "shear/tmp_shearfile2.dat"
         shear_file.save(f, [-55, 55], [30, 100, 160] , u=np.array([0.7, 1, 1.3]).T)
-        shutil.rmtree(testfilepath + "shear")
-
+        shutil.rmtree(tfp + "shear")
+        
+    def test_shear_load(self):
+        shear_file = ShearFile.load(tfp+"data/user_shear.dat")
+        np.testing.assert_array_equal(shear_file.w_positions, [30,100,160])
+        self.assertEqual(shear_file.uvw(0,65)[0],.85)
+        self.assertEqual(shear_file.uvw(-55,65)[0],.9)
+        np.testing.assert_array_equal(shear_file.uvw([0,-55],[65,65])[0],[.85,.9])
+        
+        shear_file = ShearFile.load_from_htc(tfp+"htcfiles/test.htc")
+        np.testing.assert_array_equal(shear_file.w_positions, [30,100,160])
+        np.testing.assert_array_almost_equal(shear_file.uvw([0,-55],[65,65])[0],np.array([.85,.9])*8.860807038)
 if __name__ == "__main__":
     #import sys;sys.argv = ['', 'Test.test_shearfile']
     unittest.main()