diff --git a/io/ioDymore.py b/io/ioDymore.py
new file mode 100644
index 0000000..7e75ca2
--- /dev/null
+++ b/io/ioDymore.py
@@ -0,0 +1,189 @@
+import numpy as np
+import os
+
+# ====================================================================
+#
+#
+#   READERS
+#
+#
+# ====================================================================
+
+class DymoreOutput(object):
+    def __init__(self):
+        #Different types of output
+        self.mov = None
+        self.ine = None
+        self.act = None
+        self.air = None        
+        self.frc = None
+        self.jnt = None
+        self.out = None 
+        self.freq = {}
+
+def readEigOutput():
+    '''
+    This function is to read the eigenvalues output. The returning variable is an numpy array
+    Local variables
+        :cur_dir: current path
+        :fig_dir: FIGURES path
+        :SensorEigenvalues1: First sensor eigenvalues output
+        :eig: container of eigenvalues
+    '''    
+    # cur_dir = os.getcwd()
+    # fig_dir = cur_dir+'\\FIGURES'
+    # os.chdir(fig_dir)
+    eig = np.loadtxt('SensorEigenvalues1.mdt')
+    eigvalues_last_step = eig[-1]
+    # os.chdir(cur_dir)
+    return eigvalues_last_step
+
+def readForce(fileName):
+    '''
+    fileName: <string>
+        force file name
+    '''
+    force = np.loadtxt(fileName)
+    force_last_step = force[-1]
+    return force_last_step
+	
+def makeFile(dymoreFileName):
+    """
+    Make a temporary file on the hard disk to store DYMORE-formatted data.
+
+    Parameters
+    ----------
+    dymoreFileName : <string>
+        The filename of the file to write to.
+       
+    Returns
+    -------
+    tempFile : <file object>
+        A file handle to the temporary file.
+    """
+
+    tempFile = open(dymoreFileName, 'w+')
+    return tempFile
+
+def writeDymoreMKUpdated(file, CoordType, coord, cm_x2, cm_x3, mpus, i, K, sc, gc):
+    """
+    Description.
+
+    Parameters
+    ----------
+    f : <file object>
+        The file handle that data will be written to.
+    CoordType : <string>
+        Acceptable values are: 'ETA_COORDINATE',
+                               'CURVILINEAR_COORDINATE', or
+                               'AXIAL_COORDINATE'
+    coord : <float>
+        The spanwise coordinate of this cross-section.
+        This coordinate should match the CoordType specified above.
+    cm_x2 : <float>
+        The x2-coordinate of the center of mass.
+    cm_x3 : <float>
+        The x3-coordinate of the center of mass.
+    mpus : <float>
+        The mass per unit span.
+    i1 : <float>
+        The moment of inertia about the x1-axis.
+    i2 : <float>
+        The moment of inertia about the x2-axis.
+    i3 : <float>
+        The moment of inertia about the x3-axis.
+    K : <np.array>
+        The Timoshenko stiffness matrix.
+
+    Returns
+    -------
+    <none>
+
+    Example output (written to a file)
+    ----------------------------------
+    @ETA_COORDINATE {0.00000e+00} {
+      @STIFFNESS_MATRIX { 7.6443255182E+09,   -3.5444961981E-04,   -1.5092432335E-03,    3.3599749794E+06,    2.7710007447E-01,    4.1602501550E-02,
+                                               2.8284702841E+08,   -2.8863166160E+01,   -4.5930836014E-01,   -3.3517886643E+05,    3.4162114776E-03,
+                                                                    3.5606703330E+08,   -4.0749872012E-01,    3.6079611429E-02,   -4.2577508629E+04,
+                                                                                         8.8773955810E+08,    1.8897378940E-03,    8.3869473951E-04,
+                                                                                                              4.5282893600E+10,   -5.7686739280E-02,
+                                                                                                                                   2.2625281359E+09}
+      @MASS_PER_UNIT_SPAN {7.1224712000E+02}
+      @MOMENTS_OF_INERTIA {3.9569408290E+03,
+                           3.6961203640E+03,
+                           2.6082046495E+02}
+      @CENTRE_OF_MASS_LOCATION {-1.6834618673E-17,
+                                -1.1472480873E-16}
+    }
+    """
+    f = open(file, 'w+')
+    tab = '  '
+    f.write('@BEAM_PROPERTY_DEFINITION { \n')
+    f.write(tab*2 + '@BEAM_PROPERTY_NAME {PropBeam} { \n')
+    f.write(tab*3 + '@COORDINATE_TYPE {ETA_COORDINATE} \n')
+    if CoordType == 'ETA_COORDINATE':
+        f.write(tab*2 + '@ETA_COORDINATE {' + ('%11.5e' % coord) + '} {\n')
+    elif CoordType == 'CURVILINEAR_COORDINATE':
+        # f.write(tab*2 + '@CURVILINEAR_COORDINATE {' + ('%11.5e' % coord) + '} {\n')
+        print ("***WARNING*** CURVILINEAR_COORDINATE feature is not yet supported.")
+    elif CoordType == 'AXIAL_COORDINATE':
+        f.write(tab*2 + '@AXIAL_COORDINATE {' + ('%11.5e' % coord) + '} {\n')
+    f.write(tab*3 +   '@AXIAL_STIFFNESS {' + ('%17.10e' % K[0,0])  + '} \n')
+    f.write(tab*3 +   '@BENDING_STIFFNESSES {' + ('%17.10e' % K[4,4]) + ',' + ('%20.10e' % K[5,5]) + ',' + ('%20.10e' % K[4,5])  + '} \n')   
+    f.write(tab*3 +   '@TORSIONAL_STIFFNESS {' + ('%17.10e' % K[3,3])  + '} \n')
+    f.write(tab*3 +   '@SHEARING_STIFFNESSES {' + ('%17.10e' % K[1,1]) + ',' + ('%20.10e' % K[2,2]) + ',' + ('%20.10e' % K[1,2])  + '} \n')   
+    f.write(tab*3 +   '@MASS_PER_UNIT_SPAN {' + ('%17.10e' % mpus) + '}\n')
+    f.write(tab*3 +   '@MOMENTS_OF_INERTIA {' + ('%17.10e' % i[0]) + ',\n')
+    f.write(tab*3 + ' '*21 +                    ('%17.10e' % i[1]) + ',\n')
+    f.write(tab*3 + ' '*21 +                    ('%17.10e' % i[2]) + '}\n')
+    f.write(tab*3 +   '@CENTRE_OF_MASS_LOCATION {' + ('%17.10e' % cm_x2) + ',\n')
+    f.write(tab*3 + ' '*26 +                         ('%17.10e' % cm_x3) + '}\n')
+    f.write(tab*3 +   '@SHEAR_CENTRE_LOCATION {' + ('%17.10e' % sc[1]) + ',' + ('%20.10e' % sc[2]) + '} \n')  
+    f.write(tab*3 +   '@CENTROID_LOCATION {' + ('%17.10e' % gc[1]) + ',' + ('%20.10e' % gc[2]) + '} \n')  
+    f.write(tab*2 + '}\n')
+    f.write(tab*2 + '}\n')
+    f.write('}\n')
+    f.write(tab*2 + '\n')
+    
+    return
+
+
+def calcVI(root_force, mpus, rotor_R):
+    '''
+    Compute vibration index. Followed Lim 2016, Equation 4.
+    
+    Parameters
+    
+    root_force: <array>
+        force at the root, arraged as Fx, Fy, Fz, Mx, My, Mz
+    mpus: <float>
+        mass per unit, computed from VABS
+    rotor_R:<float>
+        the radius of rotor blade
+    
+    Returns
+    
+    VI:<float>
+        vibration index
+    '''
+    
+    KF, KM = 1, 1
+    
+    W0s = mpus*rotor_R
+    W0 = W0s*4
+    
+    FxH = root_force[1]
+    FyH = root_force[2]
+    FzH = root_force[3]
+    MxH = root_force[4]
+    MyH = root_force[5]
+    
+    FH = np.sqrt((0.5*FxH)**2 + (0.67*FyH)**2 + (FzH)**2)
+    MH = np.sqrt(MxH**2 + MyH**2)
+    
+    VI = KF*FH/W0 + KM*MH/(W0*rotor_R)
+    
+    return VI
+
+if __name__ == "__main__":
+    eigvalues =readEigOutput()