diff --git a/calc_ref_state.py b/calc_ref_state.py index b0a60d5..7af441f 100644 --- a/calc_ref_state.py +++ b/calc_ref_state.py @@ -26,6 +26,7 @@ variables = {'s':{},'theta':{},'grd':[],'valid':[]} # Input 1 #################################################################################### # Physical constants +# input : None # output : r_earth, rho0, grav, gbuo # r_earth : earth radis [m] # rho0 : reference density [kg m**-3] @@ -33,10 +34,7 @@ variables = {'s':{},'theta':{},'grd':[],'valid':[]} # gbuo : # PHYSICAL CONSTATS MISSING # ############################################################################################ -r_earth = 6.4e06 -rho0 = 1027.0 -grav = 9.81 -gbuo = -grav/rho0 +# constants stored in myconstants.py # step a ##################################################################################### # Read necessary data to calculate reference state diff --git a/get_data.py b/get_data.py index a9dbb04..a61b8fc 100644 --- a/get_data.py +++ b/get_data.py @@ -1,5 +1,5 @@ -def get_data_woa2009(variables): - # get_data_woa2009(variables) ############################################################ +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 @@ -16,14 +16,19 @@ def get_data_woa2009(variables): # 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 + 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 + # 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))) # Allocation data = {} grd = [] @@ -42,35 +47,76 @@ def get_data_woa2009(variables): 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 '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'] + 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 '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'] + 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'] = np.array(DATA_d) - data[varname]['units'] = DATA_d.units - data[varname]['fill_value'] = DATA_d._FillValue + 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) if 'theta' in list_allkeys: + # compute pot. temp. for valid points + print DATA_d.dimensions + print data_perm + THETA = TEMP.transpose(data_perm) + SA = data['s']['val'] + SA = SA[data['s']['valid']] + print DATA_d.dimensions[data_perm[0]] + print DATA_d.dimensions[data_perm[1]] + print DATA_d.dimensions[data_perm[2]] + print DATA_d.dimensions[data_perm[3]] + 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 + print 'Doesnt work automatically - p_xyz' + del w3d, x3d, y3d, z3d, wvar, xvar, yvar, zvar + pr = p*0 + print 'SA.shape: ',SA.shape + print 'THETA.shape',THETA[data[varname]['valid']].shape + print 'p.shape, ',p[data[varname]['valid']].shape + dummy = input('Press enter to continue') + 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) 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 diff --git a/get_data.pyc b/get_data.pyc new file mode 100644 index 0000000..f1b2ee3 Binary files /dev/null and b/get_data.pyc differ