Calculating Vorticity, Strain in regional MITgcm output

Science
1

Hello all

I am currently using two years of hourly data from a regional MITgcm model output (configuration similar to llc4320) to estimate vorticity, strain and other kinematic quantities. I am using xmitgcm and xgcm to open the datasets stored in a server:

uv = open_mdsdataset(data_dir, grid_dir, iters='all', delta_t=150, 
                     ignore_unknown_vars=True, prefix=['diags_3D'], 
                     ref_date="2015-01-01 2:0:0", geometry='sphericalpolar',
                     llc_method='smallchunks')
grid_llc = xgcm.Grid(uv, periodic=False)

I wonder if there is any way that I could just calculate for instance, relative vorticity in the first 10 Z-levels. I have tried using the following after the two lines above:

zeta = (-grid_llc.diff(uv.UVEL.isel(Z=np.arange(0, 11)) * uv.dxC, 'Y') + 
            grid_llc.diff(uv.VVEL.isel(Z=np.arange(0, 11)) * uv.dyC, 'X'))/ds.rAz.

However, running the code yields the error

ValueError: conflicting sizes for dimension 'Z': length 10 on <this-array> and length 90 on {'time': 'time', 'XG': 'XG', 'YG': 'YG', 'Z': 'drF'}

I would like to know if there is a way to first take a slice of the model output domain and then compute zeta, instead of computing zeta for the entire domain and Z levels, which takes up to 40 hrs.

Thanks

1 Like
2

Could you show the full xarray rep for uv (output of print(uv))?

You might try slicing uv before passing it to Grid.

3

This xgcm example also might be useful for calculating vorticity efficiently

https://xgcm.readthedocs.io/en/latest/ufunc_examples.html#vorticity

4

Heres the print of (uv):

<xarray.Dataset>
Dimensions: (time: 17661, Z: 90, YC: 576, XC: 672, XG: 672, YG: 576, Zl: 90,
Zp1: 91, Zu: 90)
Coordinates: (12/33)
iter (time) int64 dask.array<chunksize=(1,), meta=np.ndarray>

  • time (time) datetime64[ns] 2015-12-27T03:00:00 … 2017-12-31T23:00:00
  • XC (XC) >f4 230.0 230.0 230.1 230.1 230.1 … 243.9 243.9 243.9 244.0
  • YC (YC) >f4 31.01 31.03 31.05 31.07 31.09 … 42.91 42.93 42.95 42.97
  • XG (XG) >f4 230.0 230.0 230.0 230.1 230.1 … 243.9 243.9 243.9 244.0
  • YG (YG) >f4 31.0 31.02 31.04 31.06 31.08 … 42.9 42.92 42.94 42.96
    … …
    maskW (Z, YC, XG) bool dask.array<chunksize=(90, 576, 672), meta=np.ndarray>
    maskS (Z, YG, XC) bool dask.array<chunksize=(90, 576, 672), meta=np.ndarray>
    maskInC (YC, XC) bool dask.array<chunksize=(576, 672), meta=np.ndarray>
    maskInS (YG, XC) bool dask.array<chunksize=(576, 672), meta=np.ndarray>
    maskInW (YC, XG) bool dask.array<chunksize=(576, 672), meta=np.ndarray>
    rhoRef (Z) >f4 dask.array<chunksize=(90,), meta=np.ndarray>
    Data variables:
    THETA (time, Z, YC, XC) float32 dask.array<chunksize=(1, 90, 576, 672), meta=np.ndarray>
    SALT (time, Z, YC, XC) float32 dask.array<chunksize=(1, 90, 576, 672), meta=np.ndarray>
    UVEL (time, Z, YC, XG) float32 dask.array<chunksize=(1, 90, 576, 672), meta=np.ndarray>
    VVEL (time, Z, YG, XC) float32 dask.array<chunksize=(1, 90, 576, 672), meta=np.ndarray>
    WVEL (time, Zl, YC, XC) float32 dask.array<chunksize=(1, 90, 576, 672), meta=np.ndarray>
    DRHODR (time, Zl, YC, XC) float32 dask.array<chunksize=(1, 90, 576, 672), meta=np.ndarray>
    Attributes:
    Conventions: CF-1.6
    title: netCDF wrapper of MITgcm MDS binary data
    source: MITgcm
    history: Created by calling `open_mdsdataset(data_dir='/data/SO2/SWO…
5

Thanks Tom. I will have a look at it.

6

I first saved UVEL and VVEL into a new xarray dataset, then sliced the new array in XC, YC, XG, YG and Z and finally passed it to Grid. Those steps really worked for me. Thanks.

3 Likes