diff --git a/weis/aeroelasticse/Turbsim_mdao/turbsim_file.py b/weis/aeroelasticse/Turbsim_mdao/turbsim_file.py
new file mode 100644
index 000000000..2756fff7b
--- /dev/null
+++ b/weis/aeroelasticse/Turbsim_mdao/turbsim_file.py
@@ -0,0 +1,319 @@
+"""Read/Write TurbSim File
+
+copied with love from the python-toolbox
+
+"""
+import pandas as pd
+import numpy as np
+import os
+import struct
+import time
+
+try:
+    from .file import File, EmptyFileError
+except:
+    EmptyFileError = type('EmptyFileError', (Exception,),{})
+    File=dict
+
+class TurbSimFile(File):
+    """ 
+    Read/write a TurbSim turbulence file (.bts). The object behaves as a dictionary.
+
+    Main keys
+    ---------
+    - 'u': velocity field, shape (3 x nt x ny x nz)
+    - 'y', 'z', 't': space and time coordinates 
+    - 'dt', 'ID', 'info'
+    - 'zTwr', 'uTwr': tower coordinates and field if present (3 x nt x nTwr)
+    - 'zHub', 'uHub': height and velocity at a reference point (usually not hub)
+
+    Main methods
+    ------------
+    - read, write, toDataFrame, keys, makePeriodic, checkPeriodic, closestPoint
+
+    Examples
+    --------
+
+        ts = TurbSimFile('Turb.bts')
+        print(ts.keys())
+        print(ts['u'].shape)  
+
+
+    """
+
+    @staticmethod
+    def defaultExtensions():
+        return ['.bts']
+
+    @staticmethod
+    def formatName():
+        return 'TurbSim binary'
+
+    def __init__(self,filename=None, **kwargs):
+        self.filename = None
+        if filename:
+            self.read(filename, **kwargs)
+
+    def read(self, filename=None, header_only=False):
+        """ read BTS file, with field: 
+                     u    (3 x nt x ny x nz)
+                     uTwr (3 x nt x nTwr)
+        """
+        if filename:
+            self.filename = filename
+        if not self.filename:
+            raise Exception('No filename provided')
+        if not os.path.isfile(self.filename):
+            raise OSError(2,'File not found:',self.filename)
+        if os.stat(self.filename).st_size == 0:
+            raise EmptyFileError('File is empty:',self.filename)
+
+        scl = np.zeros(3, np.float32); off = np.zeros(3, np.float32)
+        with open(self.filename, mode='rb') as f:            
+            # Reading header info
+            ID, nz, ny, nTwr, nt                      = struct.unpack('<h4l', f.read(2+4*4))
+            dz, dy, dt, uHub, zHub, zBottom           = struct.unpack('<6f' , f.read(6*4)  )
+            scl[0],off[0],scl[1],off[1],scl[2],off[2] = struct.unpack('<6f' , f.read(6*4))
+            nChar, = struct.unpack('<l',  f.read(4))
+            info = (f.read(nChar)).decode()
+            # Reading turbulence field
+            if not header_only: 
+                u    = np.zeros((3,nt,ny,nz))
+                uTwr = np.zeros((3,nt,nTwr))
+                # For loop on time (acts as buffer reading, and only possible way when nTwr>0)
+                for it in range(nt):
+                    Buffer = np.frombuffer(f.read(2*3*ny*nz), dtype=np.int16).astype(np.float32).reshape([3, ny, nz], order='F')
+                    u[:,it,:,:]=Buffer
+                    Buffer = np.frombuffer(f.read(2*3*nTwr), dtype=np.int16).astype(np.float32).reshape([3, nTwr], order='F')
+                    uTwr[:,it,:]=Buffer
+                u -= off[:, None, None, None]
+                u /= scl[:, None, None, None]
+                self['u']    = u
+                uTwr -= off[:, None, None]
+                uTwr /= scl[:, None, None]
+                self['uTwr'] = uTwr
+        self['info'] = info
+        self['ID']   = ID
+        self['dt']   = dt
+        self['y']    = np.arange(ny)*dy 
+        self['y']   -= np.mean(self['y']) # y always centered on 0
+        self['z']    = np.arange(nz)*dz +zBottom
+        self['t']    = np.arange(nt)*dt
+        self['zTwr'] =-np.arange(nTwr)*dz + zBottom
+        self['zHub'] = zHub
+        self['uHub'] = uHub
+
+    def write(self, filename=None):
+        """ 
+        write a BTS file, using the following keys: 'u','z','y','t','uTwr'
+                     u    (3 x nt x ny x nz)
+                     uTwr (3 x nt x nTwr)
+        """
+        if filename:
+            self.filename = filename
+        if not self.filename:
+            raise Exception('No filename provided')
+
+        nDim, nt, ny, nz = self['u'].shape
+        if 'uTwr' not in self.keys() :
+            self['uTwr']=np.zeros((3,nt,0))
+        if 'ID' not in self.keys() :
+            self['ID']=7
+
+        _, _, nTwr = self['uTwr'].shape
+        tsTwr  = self['uTwr']
+        ts     = self['u']
+        intmin = -32768
+        intrng = 65535
+        off    = np.empty((3), dtype    = np.float32)
+        scl    = np.empty((3), dtype    = np.float32)
+        info = 'Generated by TurbSimFile on {:s}.'.format(time.strftime('%d-%b-%Y at %H:%M:%S', time.localtime()))
+        # Calculate scaling, offsets and scaling data
+        out    = np.empty(ts.shape, dtype=np.int16)
+        outTwr = np.empty(tsTwr.shape, dtype=np.int16)
+        for k in range(3):
+            all_min, all_max = ts[k].min(), ts[k].max()
+            if nTwr>0:
+                all_min=min(all_min, tsTwr[k].min())
+                all_max=max(all_max, tsTwr[k].max())
+            if all_min == all_max:
+                scl[k] = 1
+            else:
+                scl[k] = intrng / (all_max-all_min)
+            off[k]    = intmin - scl[k] * all_min
+            out[k]    = (ts[k]    * scl[k] + off[k]).astype(np.int16)
+            outTwr[k] = (tsTwr[k] * scl[k] + off[k]).astype(np.int16)
+        z0 = self['z'][0]
+        dz = self['z'][1]- self['z'][0]
+        dy = self['y'][1]- self['y'][0]
+        dt = self['t'][1]- self['t'][0]
+
+        # Providing estimates of uHub and zHub even if these fields are not used
+        zHub,uHub, bHub = self.hubValues()
+
+        with open(self.filename, mode='wb') as f:            
+            f.write(struct.pack('<h4l', self['ID'], nz, ny, nTwr, nt))
+            f.write(struct.pack('<6f', dz, dy, dt, uHub, zHub, z0)) # NOTE uHub, zHub maybe not used
+            f.write(struct.pack('<6f', scl[0],off[0],scl[1],off[1],scl[2],off[2]))
+            f.write(struct.pack('<l' , len(info)))
+            f.write(info.encode())
+            for it in np.arange(nt):
+                f.write(out[:,it,:,:].tostring(order='F'))
+                f.write(outTwr[:,it,:].tostring(order='F'))
+
+    def hubValues(self, zHub=None):
+        if zHub is None:
+            try:
+                zHub=self['zHub']
+                bHub=True
+            except:
+                bHub=False
+                iz = np.argmin(np.abs(self['z']-(self['z'][0]+self['z'][-1])/2))
+                zHub = self['z'][iz]
+        else:
+            bHub=True
+        try:
+            uHub=self['uHub']
+        except:
+            iz = np.argmin(np.abs(self['z']-zHub))
+            iy = np.argmin(np.abs(self['y']-(self['y'][0]+self['y'][-1])/2))
+            uHub = np.mean(self['u'][0,:,iy,iz])
+        return zHub, uHub, bHub
+
+    def _iMid(self):
+        iy = np.argmin(np.abs(self['y']-(self['y'][0]+self['y'][-1])/2))
+        iz = np.argmin(np.abs(self['z']-(self['z'][0]+self['z'][-1])/2))
+        return iy,iz
+
+    def closestPoint(self, y, z):
+        iy = np.argmin(np.abs(self['y']-y))
+        iz = np.argmin(np.abs(self['z']-z))
+        return iy,iz
+
+    def makePeriodic(self):
+        """ Make the box periodic by mirroring it """
+        nDim, nt0, ny, nz = self['u'].shape
+        u = self['u'].copy()
+        del self['u']
+
+        nt = 2*len(self['t'])-2
+        dt = self['t'][1]- self['t'][0]
+        self['u']  = np.zeros((nDim,nt,ny,nz))
+        self['u'][:,:nt0,:,:] = u
+        self['u'][:,nt0:,:,:] = np.flip(u[:,1:-1,:,:],axis=1)
+        self['t'] = np.arange(nt)*dt
+        if 'uTwr' in self.keys():
+            _, _, nTwr = self['uTwr'].shape
+            uTwr = self['uTwr'].copy()
+            del self['uTwr']
+            # empty tower for now
+            self['uTwr'] = np.zeros((nDim,nt,nTwr))
+            self['uTwr'][:,:nt0,:] = uTwr
+            self['uTwr'][:,nt0:,:] = np.flip(uTwr[:,1:-1,:],axis=1)
+
+        self['ID']=8 # Periodic
+
+
+    def checkPeriodic(self, sigmaTol=1.5, aTol=0.5):
+        """ Check periodicity in u """
+        ic=0
+        sig  = np.std(self['u'][ic,:,:,:],axis=0)
+        mean = np.mean(self['u'][ic,:,:,:],axis=0)
+        u_first= self['u'][ic,0 ,:,:]
+        u_last = self['u'][ic,-1,:,:]
+        relSig = np.abs(u_first-u_last)/sig
+        compPeriodic = (np.max(relSig) < sigmaTol) and (np.mean(np.abs(u_first-u_last))<aTol)
+        return compPeriodic
+
+
+    def __repr__(self):
+        s='<TurbSimFile object> with keys:\n'
+        s+=' - ID {}\n'.format(self['ID'])
+        s+=' - z: [{} ... {}],  dz: {}, n: {} \n'.format(self['z'][0],self['z'][-1],self['z'][1]-self['z'][0],len(self['z']))
+        s+=' - y: [{} ... {}],  dy: {}, n: {} \n'.format(self['y'][0],self['y'][-1],self['y'][1]-self['y'][0],len(self['y']))
+        s+=' - t: [{} ... {}],  dt: {}, n: {} \n'.format(self['t'][0],self['t'][-1],self['t'][1]-self['t'][0],len(self['t']))
+        s+=' - u: ({} x {} x {} x {}) \n'.format(*(self['u'].shape))
+        ux,uy,uz=self['u'][0], self['u'][1], self['u'][2]
+        s+='    ux: min: {}, max: {}, mean: {} \n'.format(np.min(ux), np.max(ux), np.mean(ux))
+        s+='    uy: min: {}, max: {}, mean: {} \n'.format(np.min(uy), np.max(uy), np.mean(uy))
+        s+='    uz: min: {}, max: {}, mean: {} \n'.format(np.min(uz), np.max(uz), np.mean(uz))
+
+        # Mid of box, nearest neighbor
+        iy,iz = self._iMid()
+        zMid=self['z'][iz]
+        yMid=self['y'][iy]
+        uMid = np.mean(self['u'][0,:,iy,iz])
+        s+='    yMid: {} - zMid: {} - iy: {} - iz: {} - uMid: {} (nearest neighbor))\n'.format(yMid, zMid, iy, iz, uMid)
+
+#         zMid, uMid, bHub = self.hubValues()
+#         if bHub:
+#             s+='    z"Hub": {} - u"Hub": {} (NOTE: values at TurbSim "hub")\n'.format(zMid, uMid)
+
+        # Tower
+        if 'zTwr' in self.keys() and len(self['zTwr'])>0:
+            s+=' - zTwr: [{} ... {}],  dz: {}, n: {} \n'.format(self['zTwr'][0],self['zTwr'][-1],self['zTwr'][1]-self['zTwr'][0],len(self['zTwr']))
+        if 'uTwr' in self.keys() and self['uTwr'].shape[2]>0:
+            s+=' - uTwr: ({} x {} x {} ) \n'.format(*(self['uTwr'].shape))
+            ux,uy,uz=self['uTwr'][0], self['uTwr'][1], self['uTwr'][2]
+            s+='    ux: min: {}, max: {}, mean: {} \n'.format(np.min(ux), np.max(ux), np.mean(ux))
+            s+='    uy: min: {}, max: {}, mean: {} \n'.format(np.min(uy), np.max(uy), np.mean(uy))
+            s+='    uz: min: {}, max: {}, mean: {} \n'.format(np.min(uz), np.max(uz), np.mean(uz))
+            
+        return s
+
+    def toDataFrame(self):
+        dfs={}
+
+        ny = len(self['y'])
+        nz = len(self['y'])
+        # Index at mid box
+        iy,iz = self._iMid()
+
+        # Mean vertical profile
+        m = np.mean(self['u'][:,:,iy,:], axis=1)
+        s = np.std( self['u'][:,:,iy,:], axis=1)
+        ti = s/m*100
+        Cols=['z_[m]','u_[m/s]','v_[m/s]','w_[m/s]','sigma_u_[m/s]','sigma_v_[m/s]','sigma_w_[m/s]','TI_[%]']
+        data = np.column_stack((self['z'],m[0,:],m[1,:],m[2,:],s[0,:],s[1,:],s[2,:],ti[0,:]))
+        dfs['VertProfile'] = pd.DataFrame(data = data ,columns = Cols)
+
+        # Mid time series
+        u = self['u'][:,:,iy,iz]
+        Cols=['t_[s]','u_[m/s]','v_[m/s]','w_[m/s]']
+        data = np.column_stack((self['t'],u[0,:],u[1,:],u[2,:]))
+        dfs['MidLine'] = pd.DataFrame(data = data ,columns = Cols)
+
+        # Hub time series
+        #try:
+        #    zHub = self['zHub']
+        #    iz = np.argmin(np.abs(self['z']-zHub))
+        #    u = self['u'][:,:,iy,iz]
+        #    Cols=['t_[s]','u_[m/s]','v_[m/s]','w_[m/s]']
+        #    data = np.column_stack((self['t'],u[0,:],u[1,:],u[2,:]))
+        #    dfs['TSHubLine'] = pd.DataFrame(data = data ,columns = Cols)
+        #except:
+        #    pass
+        return dfs
+
+    def compute_rot_avg(self,R):
+        ''' 
+        Compute rotor average wind speed, where R is the rotor radius
+        '''
+
+        self['rot_avg'] = np.zeros((3,len(self['t'])))
+        
+        for i in range(3):
+            u_      = self['u'][i,:,:,:]
+            z_hub = self['zHub']
+            yy, zz = np.meshgrid(self['y'],self['z'])
+            rotor_ind = np.sqrt(yy**2 + (zz - z_hub)**2) < R
+
+            u_rot = []
+            for u_plane in u_:
+                u_rot.append(u_plane[rotor_ind].mean())
+
+            self['rot_avg'][i,:] = u_rot
+
+
+if __name__=='__main__':
+    ts = TurbSimFile('../_tests/TurbSim.bts')
diff --git a/weis/aeroelasticse/Turbsim_mdao/turbsim_writer.py b/weis/aeroelasticse/Turbsim_mdao/turbsim_writer.py
index 20a041e90..d5253a534 100644
--- a/weis/aeroelasticse/Turbsim_mdao/turbsim_writer.py
+++ b/weis/aeroelasticse/Turbsim_mdao/turbsim_writer.py
@@ -1,5 +1,4 @@
 from weis.aeroelasticse.Turbsim_mdao.turbsim_vartrees import turbsiminputs
-from weis.aeroelasticse.Turbsim_mdao.turbulence_spectrum import turb_specs
 from weis.aeroelasticse.Turbsim_mdao.wind_profile_writer import write_wind
 import os
 import numpy as np
@@ -14,14 +13,7 @@ def __init__(self):
  
          # Turbulence file parameters
          self.wind_speed = 8.
-         self.L_u = 2.54e+02
-         self.L_v=1.635e+02
-         self.L_w=4.7e+01
-         self.sigma_u=1.325
-         self.sigma_v=0.9
-         self.sigma_w=0.7625
          self.turbulence_file_name = 'tsim_user_turbulence_default.inp'
-         self.turbulence_template_file = 'TurbsimInputFiles/turbulence_user.inp'
          
          # profile file parameters
          self.profile_template = 'TurbsimInputFiles/shear.profile'
@@ -31,7 +23,7 @@ def __init__(self):
 
          self.run_dir = 'run%d'%np.random.uniform(0,1e10)
 
-    def execute(self, write_specs=False, write_profile=False):
+    def execute(self, write_profile=False):
          if not os.path.exists(self.run_dir): 
             try:
                sleep(random.uniform(0, 1))
@@ -40,12 +32,6 @@ def execute(self, write_specs=False, write_profile=False):
             except:
                pass
 
-         # Write turbulence file
-         if write_specs:
-            self.turbsim_vt.metboundconds.UserFile = os.sep.join([self.run_dir, self.turbulence_file_name])
-            turb_specs(V_ref=float(self.wind_speed), L_u=float(self.L_u), L_v=float(self.L_v), L_w=float(self.L_w), sigma_u=float(self.sigma_u), sigma_v=float(self.sigma_v), sigma_w=float(self.sigma_w), filename=self.turbsim_vt.metboundconds.UserFile, template_file=self.turbulence_template_file)
-            self.turbsim_vt.metboundconds.UserFile = os.sep.join(['..', self.run_dir, self.turbulence_file_name])
-
          # Write profile file
          if write_profile:
             self.turbsim_vt.metboundconds.ProfileFile = os.sep.join([self.run_dir, self.turbsim_vt.metboundconds.ProfileFile])
diff --git a/weis/aeroelasticse/Turbsim_mdao/turbulence_spectrum.py b/weis/aeroelasticse/Turbsim_mdao/turbulence_spectrum.py
index ca4e3ab02..b10b345b2 100644
--- a/weis/aeroelasticse/Turbsim_mdao/turbulence_spectrum.py
+++ b/weis/aeroelasticse/Turbsim_mdao/turbulence_spectrum.py
@@ -1,37 +1,133 @@
 import numpy as np
-import os
-# import pandas
-from collections import OrderedDict
-import weis.aeroelasticse
-def turb_specs(V_ref, L_u, L_v, L_w, sigma_u, sigma_v, sigma_w, template_file, filename):
-    
-    f=(np.array([np.arange(0.0015873015873015873015873015873, 20.00001, 0.0015873015873015873015873015873)])).T
-    
-    a=int(len(f))
-    
-    U=np.zeros((a,1),dtype=float)
-    V=np.zeros((a,1),dtype=float)
-    W=np.zeros((a,1),dtype=float)
-    
-    for i in range(0,a):
-        U[i,0]= (4*L_u/V_ref)*sigma_u**2/((1+6*f[i,0]*L_u/V_ref)**(5./3.))
-        V[i,0]= (4*L_v/V_ref)*sigma_v**2/((1+6*f[i,0]*L_v/V_ref)**(5./3.))
-        W[i,0]= (4*L_w/V_ref)*sigma_w**2/((1+6*f[i,0]*L_w/V_ref)**(5./3.))
-        
-    df=pandas.DataFrame({'Frequency (Hz)':f[:,0],'u-component PSD (m^2/s)': U[:,0],'v-component PSD (m^2/s)': V[:,0],'w-component PSD (m^2/s)':W[:,0]})
-    with open(template_file, 'r') as f:
-        get_all=f.readlines() 
+import matplotlib.pyplot as plt
+from weis.aeroelasticse.pyIECWind import pyIECWind_extreme
+from scipy.special import modstruve, iv
+
+from ROSCO_toolbox.ofTools.util import spectral
+from weis.aeroelasticse.Turbsim_mdao.turbsim_file import TurbSimFile
+
+
+def IECKaimal(f, V_ref, HH, Class, Categ, TurbMod, R):
     
+    ###### Initialize IEC Wind parameters #######
+    iec_wind = pyIECWind_extreme()
+    iec_wind.z_hub = HH
+    iec_wind.Turbine_Class = Class
+    iec_wind.Turbulence_Class = Categ
+    iec_wind.setup()
+
+    # Compute wind turbulence standard deviation (invariant with height)
+    if TurbMod == 'NTM':
+        sigma_1 = iec_wind.NTM(V_ref)
+    elif TurbMod == 'ETM':
+        sigma_1 = iec_wind.ETM(V_ref)
+    elif TurbMod == 'EWM':
+        sigma_1 = iec_wind.EWM(V_ref)
+    else:
+        raise Exception("Wind model must be either NTM, ETM, or EWM. While you wrote " + TurbMod)
+
+    # Compute turbulence scale parameter Annex C3 of IEC 61400-1-2019
+    # Longitudinal
+    if HH <= 60:
+        L_1 = .7 * HH
+    else:
+        L_1 = 42.
+    sigma_u = sigma_1
+    L_u = 8.1 * L_1
+    # Lateral
+    sigma_v =  0.8 * sigma_1
+    L_v = 2.7 * L_1 
+    # Upward
+    sigma_w =  0.5 * sigma_1
+    L_w = 0.66 * L_1 
+
+    U = (4*L_u/V_ref)*sigma_u**2/((1+6*f*L_u/V_ref)**(5./3.))
+    V = (4*L_v/V_ref)*sigma_v**2/((1+6*f*L_v/V_ref)**(5./3.))
+    W = (4*L_w/V_ref)*sigma_w**2/((1+6*f*L_w/V_ref)**(5./3.))
+
+    kappa = 12 * np.sqrt((f/V_ref)**2 + (0.12 / L_u)**2)
+
+    Rot = (2*U / (R * kappa)**3) * \
+        (modstruve(1,2*R*kappa) - iv(1,2*R*kappa) - 2/np.pi + \
+            R*kappa * (-2 * modstruve(-2,2*R*kappa) + 2 * iv(2,2*R*kappa) + 1) )
+
+    # set NaNs to 0
+    Rot[np.isnan(Rot)] = 0
     
+    # Formulas from Section 6.3 of IEC 61400-1-2019
+    # S_1_f = 0.05 * sigma_1**2. * (L_1 / V_hub) ** (-2./3.) * f **(-5./3)
+    # S_2_f = S_3_f = 4. / 3. * S_1_f
+    # sigma_k = np.sqrt(np.trapz(S_1_f, f))
+    # print(sigma_k)
+    # print(sigma_u)
+
+    return U, V, W, Rot
+
+if __name__=="__main__":
+
+    # Average wind speed at hub height
+    V_hub = 6.
+    # Wind turbine hub height
+    HH = 150.
+    # IEC Turbine Wind Speed Class, can be I, II, or III
+    Class = 'I'
+    # IEC Turbine Wind Turbulence Category, can be A, B, or C
+    Categ = 'B'
+    # Wind turbulence model, it can be NTM = normal turbulence, ETM = extreme turbulence, or EWM = extreme wind
+    TurbMod = 'NTM'
+    # Rotor radius, used to calculate rotor average wind spectrum
+    R = 120.
+    # Frequency range
+    f=np.arange(0.0015873015873015873015873015873, 20.00001, 0.0015873015873015873015873015873)
+    f=np.logspace(-2,1)
+
+    U, V, W, S_r = IECKaimal(f, V_hub, HH, Class, Categ, TurbMod, R)
+
+
+
+
+    # load turbsim file and compute spectra of rot avg wind speed
+    # ts_file = TurbSimFile('/Users/dzalkind/Tools/WEIS-2/wind/IEA-15MW/level3_NTM_U12.000000_Seed600.0.bts')
+    ts_filename   = '/Users/dzalkind/Tools/WEIS-2/wind/IEA-15MW/level3_NTM_U6.000000_Seed600.0.bts'
+    ts_file = TurbSimFile(ts_filename)
+    ts_file.compute_rot_avg(120)
+
+    dist = {}
+    dist['Time'] = ts_file['t']
+    dist['Wind'] = ts_file['rot_avg'][0]
+    dist['HH'] = ts_file['u'][0,:,12,12]
+
+
+    # y,fq= op.plot_spectral([dist],[('Wind',0)],averaging='Welch',averaging_window='Hamming')
+    fq, y, _ = spectral.fft_wrap(
+                    dist['Time'], dist['Wind'], averaging='Welch', averaging_window='Hamming', output_type='psd')
+
+    plt.plot(fq, (y), '-', label = 'U_TS')
+
+    plt.plot(f, U, '-', label = 'U')
+    plt.plot(f, V, '--', label = 'V')
+    plt.plot(f, W, ':', label = 'W')
+    plt.plot(f, S_r, label = 'S_rot')
+    plt.yscale('log')
+    plt.xscale('log')
+
+    plt.xlim([1e-2,10])
+    plt.grid('True')
+
+    plt.xlabel('Freq. (Hz)')
+    plt.ylabel('PSD')
+
+
+    # How much energy are we capturing with f < 0.3
+    if False:
+        f = f[~np.isnan(S_r)]
+        s_r = np.sqrt(S_r[~np.isnan(S_r)])
+
+        tot_eng = np.trapz(s_r,f)
+
+        eng = np.array([np.trapz(s_r[:i],f[:i]) for i in range(len(f))]) / tot_eng
         
-    with open(filename,'w') as f:
-        for i,line in enumerate(get_all):
-                if i < 11:  ## STARTS THE NUMBERING FROM 1 (by default it begins with 0)
-                        f.writelines(line)                             ## OVERWRITES line:2
-                else: 
-                        f.write(df.to_string(index=False,header=False,col_space=15))
-                        break
-                
-                                                                                                               
-    
+        plt.plot(f,eng * 1e2,label='Perc. Energy')
 
+    plt.legend()
+    plt.show() 
\ No newline at end of file
diff --git a/weis/aeroelasticse/test/test_IECWind.py b/weis/aeroelasticse/test/test_IECWind.py
index 7987d1f2b..bbe7a7bd2 100644
--- a/weis/aeroelasticse/test/test_IECWind.py
+++ b/weis/aeroelasticse/test/test_IECWind.py
@@ -1,8 +1,7 @@
 import unittest
-
-from openmdao.utils.assert_utils import assert_near_equal
 from weis.aeroelasticse.pyIECWind import pyIECWind_extreme
-
+from weis.aeroelasticse.Turbsim_mdao.turbulence_spectrum import IECKaimal
+import numpy as np
 
 class TestIECWind(unittest.TestCase):
     def setUp(self):
@@ -11,20 +10,55 @@ def setUp(self):
 
     def test_NTM(self):
         sigma_1 = self.extreme_wind.NTM(10.0)
-        assert_near_equal(sigma_1, 1.834)
+        np.testing.assert_almost_equal(sigma_1, 1.834)
 
     def test_ETM(self):
         sigma_1 = self.extreme_wind.ETM(10.0)
-        assert_near_equal(sigma_1, 2.96128)
+        np.testing.assert_almost_equal(sigma_1, 2.96128)
 
     def test_EWM(self):
         sigma_1, V_e50, V_e1, V_50, V_1 = self.extreme_wind.EWM(10.0)
-        assert_near_equal(sigma_1, 1.1)
-        assert_near_equal(V_e50, 70.0)
-        assert_near_equal(V_e1, 56.0)
-        assert_near_equal(V_50, 50.0)
-        assert_near_equal(V_1, 40.0)
+        np.testing.assert_almost_equal(sigma_1, 1.1)
+        np.testing.assert_almost_equal(V_e50, 70.0)
+        np.testing.assert_almost_equal(V_e1, 56.0)
+        np.testing.assert_almost_equal(V_50, 50.0)
+        np.testing.assert_almost_equal(V_1, 40.0)
+
+    def testIECKaimal(self):
+        # Average wind speed at hub height
+        V_hub = 10.
+        # Wind turbine hub height
+        HH = 100.
+        # IEC Turbine Wind Speed Class, can be I, II, or III
+        Class = 'I'
+        # IEC Turbine Wind Turbulence Category, can be A, B, or C
+        Categ = 'A'
+        # Wind turbulence model, it can be NTM = normal turbulence, ETM = extreme turbulence, or EWM = extreme wind
+        TurbMod = 'NTM'
+        # Rotor radius, used to calculate rotor average wind spectrum
+        R = 60.
+        # Frequency range
+        f=np.linspace(0.001, 20., 10)
+
 
+        U, V, W, Rot = IECKaimal(f, V_hub, HH, Class, Categ, TurbMod, R)
+        np.testing.assert_array_almost_equal(U , np.array([4.38659198e+02, 2.22282686e-02, 7.01694129e-03, 
+                                                           3.57258290e-03, 2.21264170e-03, 1.52577475e-03, 
+                                                           1.12612267e-03, 8.71070791e-04, 6.97323831e-04,
+                                                           5.73069089e-04]), 5)
+        np.testing.assert_array_almost_equal(V , np.array([1.14285126e+02, 2.93758372e-02, 9.30714989e-03,
+                                                           4.74439644e-03, 2.94018708e-03, 2.02821285e-03, 
+                                                           1.49732131e-03, 1.15840022e-03, 9.27463030e-04, 
+                                                           7.62277955e-04]), 5)
+        np.testing.assert_array_almost_equal(W , np.array([1.18477575e+01, 2.83846534e-02, 9.14368382e-03,
+                                                           4.68723385e-03, 2.91293916e-03, 2.01281628e-03,
+                                                           1.48763283e-03, 1.15183382e-03, 9.22764383e-04,
+                                                           7.58773687e-04]), 5)
+        np.testing.assert_array_almost_equal(Rot , np.array([3.47567168e+02, 1.72828935e-06, 1.36729085e-07,
+                                                            0,           0,            0,    
+                                                            0,           0,            0,
+                                                            0]), 5)
+                                                           
 
 if __name__ == "__main__":
     unittest.main()