Merge pull request #38 from andrewdelman/geos_bal_tutorial

Geostrophic balance tutorial

Author: andrewdelman

ECCO-ACCESS/Downloading_ECCO_datasets_from_PODAAC/Tutorial_Python3_Jupyter_Notebook_Downloading_ECCO_Datasets_from_PODAAC.ipynb

@@ -925,7 +925,7 @@
     "\n",
     "[**ecco_download** module](../../ECCO-ACCESS/Downloading_ECCO_datasets_from_PODAAC/ecco_download.py)\n",
     "\n",
-    "You can save this file either in the same directory where you store the tutorial notebooks, or a different directory that you then add to your path using sys.path.append. The syntax for using this module will be discussed in the next tutorial (Geostrophic Balance)."
+    "You can save this file either in the same directory where you store the tutorial notebooks, or a different directory that you then add to your path using sys.path.append. You can see the syntax of how the module is used in upcoming tutorials."
    ]
   }
  ],

Intro_to_PO_Tutorials/Geostrophic_balance.ipynb

@@ -0,0 +1 @@
+../ECCO-ACCESS/Downloading_ECCO_datasets_from_PODAAC/ecco_download.py

Intro_to_PO_Tutorials/ecco_download.py

@@ -0,0 +1,78 @@
+def geos_vel_compute(dens_press_filename,grid_filename="~/Downloads/ECCO_V4r4_PODAAC/ECCO_L4_GEOMETRY_LLC0090GRID_V4R4/GRID_GEOMETRY_ECCO_V4r4_native_llc0090.nc",fc_filename="llc_13tile_fc.txt"):
+    """
+    This routine computes geostrophic velocities from an input netCDF file containing ECCO v4r4 density and pressure anomalies on the native llc90 grid.
+    
+    Parameters
+    ----------
+    dens_press_filename: the name (including path if not in current directory) of the netCDF file containing the density and pressure anomalies
+    
+    grid_filename: the name (including path if not in current directory) of the netCDF file containing the latitudes ('YC') and grid parameters ('dxC','dyC') needed to compute horizontal derivatives
+    
+    fc_filename: text file containing a Python dictionary specifying the tile/face connections in the ECCO llc90 grid
+    """
+    
+    import numpy as np
+    import xarray as xr
+    import json
+    import xgcm
+    
+    # load file into workspace
+    ds_denspress = xr.open_dataset(dens_press_filename, data_vars='minimal',\
+                                     coords='minimal', compat='override')
+    
+    densanom = ds_denspress.RHOAnoma
+    
+    rhoConst = 1029.
+    dens = rhoConst + densanom
+    
+    pressanom = ds_denspress.PHIHYDcR
+    
+    
+    # load grid parameters file
+    ds_grid = xr.open_dataset(grid_filename)
+    
+    
+    # load face_connections dictionary
+    
+    # read in as string
+    with open(fc_filename) as fc:
+        data = fc.read()
+    # convert string to dictionary
+    face_connections = json.loads(data)
+    
+    
+    # create xgcm Grid object
+    grid = xgcm.Grid(ds_grid,periodic=False,face_connections=face_connections)
+    
+    # compute derivatives of pressure in x and y
+    d_press_dx = (grid.diff(rhoConst*pressanom,axis="X",boundary='extend'))/ds_grid.dxC
+    d_press_dy = (grid.diff(rhoConst*pressanom,axis="Y",boundary='extend'))/ds_grid.dyC
+    
+    # interpolate (vector) gradient values to center of grid cells
+    press_grads_interp = grid.interp_2d_vector({'X':d_press_dx,'Y':d_press_dy},boundary='extend')
+    dp_dx = press_grads_interp['X']
+    dp_dy = press_grads_interp['Y']
+    
+    # compute f from latitude of grid cell centers
+    lat = ds_grid.YC
+    Omega = (2*np.pi)/86164
+    lat_rad = (np.pi/180)*lat    # convert latitude from degrees to radians
+    f = 2*Omega*np.sin(lat_rad)
+    
+    # compute geostrophic velocities
+    v_g = dp_dx/(f*dens)
+    u_g = -dp_dy/(f*dens)
+    
+    # assign attributes to DataArrays (names) and units
+    u_g.name = 'u_g'
+    u_g.attrs.update({'long_name': 'Geostrophic velocity in model-x direction',\
+                      'units': 'm s-1'})
+    v_g.name = 'v_g'
+    v_g.attrs.update({'long_name': 'Geostrophic velocity in model-y direction',\
+                      'units': 'm s-1'})
+    
+    # create xarray Dataset containing geostrophic velocities
+    ds_geos_vel = u_g.to_dataset(name='u_g',promote_attrs=True)
+    ds_geos_vel['v_g'] = v_g
+    
+    return ds_geos_vel

Intro_to_PO_Tutorials/ecco_download.py

@@ -0,0 +1 @@
+{'tile': {0: {'X': ((12, 'Y', False), (3, 'X', False)), 'Y': (None, (1, 'Y', False))}, 1: {'X': ((11, 'Y', False), (4, 'X', False)), 'Y': ((0, 'Y', False), (2, 'Y', False))}, 2: {'X': ((10, 'Y', False), (5, 'X', False)), 'Y': ((1, 'Y', False), (6, 'X', False))}, 3: {'X': ((0, 'X', False), (9, 'Y', False)), 'Y': (None, (4, 'Y', False))}, 4: {'X': ((1, 'X', False), (8, 'Y', False)), 'Y': ((3, 'Y', False), (5, 'Y', False))}, 5: {'X': ((2, 'X', False), (7, 'Y', False)), 'Y': ((4, 'Y', False), (6, 'Y', False))}, 6: {'X': ((2, 'Y', False), (7, 'X', False)), 'Y': ((5, 'Y', False), (10, 'X', False))}, 7: {'X': ((6, 'X', False), (8, 'X', False)), 'Y': ((5, 'X', False), (10, 'Y', False))}, 8: {'X': ((7, 'X', False), (9, 'X', False)), 'Y': ((4, 'X', False), (11, 'Y', False))}, 9: {'X': ((8, 'X', False), None), 'Y': ((3, 'X', False), (12, 'Y', False))}, 10: {'X': ((6, 'Y', False), (11, 'X', False)), 'Y': ((7, 'Y', False), (2, 'X', False))}, 11: {'X': ((10, 'X', False), (12, 'X', False)), 'Y': ((8, 'Y', False), (1, 'X', False))}, 12: {'X': ((11, 'X', False), None), 'Y': ((9, 'Y', False), (0, 'X', False))}}}

Intro_to_PO_Tutorials/geos_vel.py

@@ -0,0 +1 @@
+../Intro_to_PO_Tutorials/Geostrophic_balance.ipynb

Intro_to_PO_Tutorials/llc_13tile_fc.txt

@@ -70,6 +70,12 @@ The `ecco_v4_py`_ package used in this tutorial was inspired by the `xmitgcm`_ p
 
    ECCO_v4_Example_calculations_with_scalar_quantities.ipynb
 
+.. toctree::
+   :maxdepth: 2
+   :caption: Intro to PO Tutorials
+   
+   Geostrophic_balance.ipynb
+
 .. toctree::
    :maxdepth: 2
    :caption: More Advanced Calculations