def get_data_woa2009(variables,grid_order=('lon','lat','z','time',)):
    # get_data_woa2009(variables,grid_order) #################################################
    # 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 '  Load WOA2009 data'
    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  # class containing grid information
    import myconstants as my_const   # my own defined constants
    from mycalc_ocean import calc_theta_from_temp

    # Allocation and Initialization %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    # define grid_order (used for data permutation) 
    find_index = lambda searchlist, elem: [[i for i, x in enumerate(searchlist) if x == e] for e in elem]
    grid_order = np.squeeze(np.array(find_index(('lon','lat','z','time',),grid_order)))
    list_allkeys = variables.keys()
    # Allocation
    data = {}
    grd  = []
    
    # Read grid information (for all variables identical) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    print '  -> Read grid information'
    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)
    # define data permutation to match input dimensions
    nc_fid.variables['t_an'].dimensions
    
    # read salinity data %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if 's' in list_allkeys:
        print '  -> Read salinity data'
        fn        = dir_data + 'salinity_annual_1deg.nc'
        nc_fid    = Dataset(fn,'r')
        # read salinity data
        DATA_d    = nc_fid.variables['s_an']
        MASK      = np.array(DATA_d)!=DATA_d._FillValue
        data_perm = np.squeeze(np.asarray(find_index(DATA_d.dimensions,[grid_names[i] for i in grid_order])))
        varname   = 's'
        data[varname]                  = {}
        data[varname]['val']           = np.array(DATA_d).transpose(data_perm)
        data[varname]['units']         = DATA_d.units
        data[varname]['fill_value']    = DATA_d._FillValue
        data[varname]['standard_name'] = DATA_d.standard_name
        data[varname]['valid']         = MASK.transpose(data_perm)

    # Read temperature data %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if 'temp' or 'theta' in list_allkeys:
        print '  -> Read temperature data'
        fn       = dir_data + 'temperature_annual_1deg.nc'
        nc_fid   = Dataset(fn,'r')
        # read temperature data
        DATA_d    = nc_fid.variables['t_an']
        data_perm = np.squeeze(np.asarray(find_index(DATA_d.dimensions,[grid_names[i] for i in grid_order])))
        TEMP      = np.array(DATA_d)
        MASK      = TEMP!=DATA_d._FillValue
        if 'temp' in list_allkeys:
            varname = 'temp'
            data[varname]                  = {}
            data[varname]['val']           = TEMP.transpose(data_perm)
            data[varname]['units']         = DATA_d.units
            data[varname]['fill_value']    = DATA_d._FillValue
            data[varname]['standard_name'] = DATA_d.standard_name
            data[varname]['valid']         = MASK.transpose(data_perm)

        # compute potential temperature for valid points -------------------------------------
        if 'theta' in list_allkeys:
            print '  -> Compute potential temperature from temperature'
            THETA       = np.copy(TEMP).transpose(data_perm)
            SA          = data['s']['val']
            SA          = SA[data['s']['valid']]
            wvar        = nc_fid.variables[DATA_d.dimensions[data_perm[0]]]
            xvar        = nc_fid.variables[DATA_d.dimensions[data_perm[1]]]
            yvar        = nc_fid.variables[DATA_d.dimensions[data_perm[2]]]
            zvar        = nc_fid.variables[DATA_d.dimensions[data_perm[3]]]
            w3d, x3d, y3d, z3d = np.meshgrid(xvar,wvar,yvar,zvar)
            p           = y3d*(my_const.rho0*my_const.grav/1e4) # in dbar
            del w3d, x3d, y3d, z3d, wvar, xvar, yvar, zvar
            pr          = p*0.0
            print '  *** why do I set p_ref=0? ***'
            MASK_mesh   = data[varname]['valid']
            THETA[MASK_mesh] = calc_theta_from_temp(SA,THETA[MASK_mesh],\
                p[MASK_mesh],pr[MASK_mesh])
            varname     = 'theta'
            data[varname]                  = {}
            data[varname]['val']           = THETA
            data[varname]['units']         = DATA_d.units
            data[varname]['fill_value']    = DATA_d._FillValue
            data[varname]['standard_name'] = 'sea_water_potential_temperature'
            data[varname]['valid']         = MASK.transpose(data_perm)

    # return data %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if 'grd' in list_allkeys:
        return grd, data
    else:
        return data