Update of get_data and calc_ref_state: running version to read data

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
+# constants stored in myconstants.py
-rho0    = 1027.0
-grav    = 9.81
-gbuo    = -grav/rho0
 
 # step a #####################################################################################
 # Read necessary data to calculate reference state

get_data.py

@@ -1,5 +1,5 @@
-def get_data_woa2009(variables):
+def get_data_woa2009(variables,grid_order=('lon','lat','z','time',)):
-    # get_data_woa2009(variables) ############################################################
+    # 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 netCDF4 import Dataset      # to read netcdf files
-    from mydata_classes import Grid
+    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]                  = {}
-            data[varname]['val']        = np.array(DATA_d)
+            data[varname]['val']           = TEMP.transpose(data_perm)
-            data[varname]['units']      = DATA_d.units
+            data[varname]['units']         = DATA_d.units
-            data[varname]['fill_value'] = DATA_d._FillValue
+            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