diff --git a/calc_ref_state.py b/calc_ref_state.py
index 475064f..974eba7 100644
--- a/calc_ref_state.py
+++ b/calc_ref_state.py
@@ -8,16 +8,20 @@
 # Load necessary modules #####################################################################
 import numpy as np      # 
 import time             # to measure elapsed time
+from get_data import get_data_woa2009 as read_data
 # ############################################################################################
 
 # Input 0 ####################################################################################
 # Manual defined inputs
 # input         : None 
-# output        : dir_plot
-#     dir_plot  : directory, where plots will be saved 
+# output        : dir_plot, variables
+#     dir_plot  : directory, where plots will be saved
+#     variables : variables necessary for calculations 
 # MANUAL DEFINED OUTPUTS MISSING
 # ############################################################################################
 print 'MANUAL DEFINED OUTPUTS MISSING'
+dir_plot  = '/glusterfs/inspect/users/xg911182/Code/Python/plots_test/'
+variables = {'s':{},'theta':{},'grd':[],'valid':[]}
 
 # Input 1 ####################################################################################
 # Physical constants
@@ -41,7 +45,8 @@ print 'PHYSICAL CONSTATS MISSING'
 beg_time = time.time()
 print ''
 print ' ----------------------------------------------'
-print 'Insert code to read data'
+print ' Read data'
+grd, data = read_data(variables)
 print ' ----------------------------------------------'
 end_time = time.time()
 print 'Elapsed time to read data:         '+ '{:.2f}'.format(end_time-beg_time)+'s'
diff --git a/get_data.py b/get_data.py
index 30f72c5..a9dbb04 100644
--- a/get_data.py
+++ b/get_data.py
@@ -17,10 +17,62 @@ def get_data_woa2009(variables):
     #                    Ed. NOAA Atlas NESDIS 69, 184 pp.
     #   -> PDF         : https://www.nodc.noaa.gov/OC5/indpub.html#woa09
     # #######################################################################################
-    dir_data = '/glusterfs/inspect/users/xg911182/data/WOA2009/'
+    print 'Use of get_data_woa2009 in get_data.py'
+    dir_data   = '/glusterfs/inspect/users/xg911182/data/WOA2009/'
+    grid_names = ('lon','lat','depth','time',)
+    # import modules
+    import numpy as np
+    from netCDF4 import Dataset # to read netcdf files
+    from mydata_classes import Grid
     # Allocation
     data = {}
+    grd  = []
+    list_allkeys = variables.keys()
+    # Read grid information (for all variables identical)
+    print 'Check if all WOA data is using the same grid'
+    if 'grd' in list_allkeys:
+        fn       = dir_data + 'temperature_annual_1deg.nc'
+        nc_fid   = Dataset(fn,'r')
+        LON      = np.array(nc_fid.variables[grid_names[0]])
+        LAT      = np.array(nc_fid.variables[grid_names[1]])
+        Z        = np.array(nc_fid.variables[grid_names[2]])
+        TIME     = nc_fid.variables[grid_names[3]]
+        Ulon     = nc_fid.variables[grid_names[0]].units
+        Ulat     = nc_fid.variables[grid_names[1]].units
+        Uz       = nc_fid.variables[grid_names[2]].units
+        Ut       = nc_fid.variables[grid_names[3]].units
+        grd      = Grid(LON,LAT,Z,TIME,len(LON),len(LAT),len(Z),len(TIME), Ulon, Ulat, Uz, Ut)
     # Read temperature data
+    if 'temp' or 'theta' in list_allkeys:
+        print 'Load temperature data'
+        fn       = dir_data + 'temperature_annual_1deg.nc'
+        nc_fid   = Dataset(fn,'r')
+        # read temperature data
+        DATA_d   = nc_fid.variables['t_an']
+        if 'temp' in list_allkeys:
+            varname = 'temp'
+            data[varname]               = {}
+            data[varname]['val']        = np.array(DATA_d)
+            data[varname]['units']      = DATA_d.units
+            data[varname]['fill_value'] = DATA_d._FillValue
+        if 'theta' in list_allkeys:
+            print 'MISSING: Calculate pot. temp. from temperature data'
     # read salinity data
+    if 's' in list_allkeys:
+        print 'Load salinity data'
+        fn      = dir_data + 'salinity_annual_1deg.nc'
+        nc_fid  = Dataset(fn,'r')
+        # read salinity data
+        DATA_d  = nc_fid.variables['s_an']
+        varname = 's'
+        data[varname]               = {}
+        data[varname]['val']        = np.array(DATA_d)
+        data[varname]['units']      = DATA_d.units
+        data[varname]['fill_value'] = DATA_d._FillValue
+
+
     # return data
-    return data 
+    if 'grd' in list_allkeys:
+        return grd, data
+    else:
+        return data