def get_data_woa2009(variables):
    # get_data_woa2009(variables) ############################################################
    # written by    : Gabriel Wolf, g.a.wolf@reading.ac.uk
    #                 adapted from get_woa2009_data.m of Remi Tailleux
    # last modified : 13.09.2018
    # Content/Description ####################################################################
    # Loading data from the 2009 version of the Levitus World Ocean Atlas
    # References:
    # - for temperature: Locarnini, R. A., A. V. Mishonov, J. I. Antonov, T. P. Boyer, H. E. 
    #                    Garcia, O. K. Baranova, M. M. Zweng, and D. R. Johnson, 2010: World 
    #                    Ocean Atlas 2009, Volume 1: Temperature. S. Levitus, Ed. NOAA Atlas
    #                    NESDIS 68, 184 pp.
    #   -> PDF         : https://www.nodc.noaa.gov/OC5/indpub.html#woa09
    # - for salinity   : Antonov, J. I., D. Seidov, T. P. Boyer, R. A. Locarnini, A. V. 
    #                    Mishonov, H. E. Garcia, O. K. Baranova, M. M. Zweng, and D. R. 
    #                    Johnson, 2010: World Ocean Atlas 2009, Volume 2: Salinity. S. Levitus,
    #                    Ed. NOAA Atlas NESDIS 69, 184 pp.
    #   -> PDF         : https://www.nodc.noaa.gov/OC5/indpub.html#woa09
    # #######################################################################################
    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
    if 'grd' in list_allkeys:
        return grd, data
    else:
        return data