Reduced Shear Evolution equation#

\begin{equation} \frac{D}{Dt} Sh^2_{red} = u_z^2v_y +v_z^2u_x - (u_y + v_x)u_zv_z + u_z F_z^x + v_z F_z^y + \left(\frac{1}{2Ri_c - 1}\right) (u_zb_x + v_zb_y) + \frac{1}{2Ri_c} w_zb_z + \frac{1}{2Ri_c}∂_{zz}Kb_z \end{equation}

# %load_ext %autoreload
%load_ext watermark
%autoreload 1

%matplotlib inline
%config InlineBackend.figure_format="retina"

import cf_xarray
import dask
import dcpy
import distributed
import holoviews as hv
import hvplot.xarray
import matplotlib as mpl
import matplotlib.dates as mdates
import matplotlib.pyplot as plt
import matplotlib.units as munits
import numpy as np
import pandas as pd
import panel.widgets as pnw
import seawater as sw
import xarray as xr
import xgcm
from holoviews import opts

%aimport pump


mpl.rcParams["figure.dpi"] = 140
mpl.rcParams["savefig.dpi"] = 300
mpl.rcParams["savefig.bbox"] = "tight"

munits.registry[np.datetime64] = mdates.ConciseDateConverter()

xr.set_options(keep_attrs=True)


hv.extension("bokeh")
# hv.opts.defaults(opts.Image(fontscale=1.5), opts.Curve(fontscale=1.5))


%watermark -iv


xr.DataArray([1.0])

panel      : 0.12.0
seawater   : 3.3.4
dask       : 2021.7.2
cf_xarray  : 0.6.0
dcpy       : 0.1
xgcm       : 0.5.1.dev138+g3b80993
distributed: 2021.7.2
hvplot     : 0.7.3
xarray     : 0.19.0
pump       : 0.1
matplotlib : 3.4.2
pandas     : 1.3.1
holoviews  : 1.14.5
numpy      : 1.21.1

<xarray.DataArray (dim_0: 1)>
array([1.])
Dimensions without coordinates: dim_0
import dask_jobqueue

if "client" in locals():
    client.close()
# if "cluster" in locals():
#    cluster.close()

env = {"OMP_NUM_THREADS": "1", "NUMBA_NUM_THREADS": "1", "OPENBLAS_NUM_THREADS": "1"}

# cluster = distributed.LocalCluster(
#    n_workers=8,
#    threads_per_worker=1,
#    env=env
# )

if "cluster" in locals():
    del cluster

# cluster = ncar_jobqueue.NCARCluster(
#    account="NCGD0011",
#    scheduler_options=dict(dashboard_address=":9797"),
#    env_extra=env,
# )
cluster = dask_jobqueue.PBSCluster(
    cores=1,  # The number of cores you want
    memory="23GB",  # Amount of memory
    processes=1,  # How many processes
    queue="casper",  # The type of queue to utilize (/glade/u/apps/dav/opt/usr/bin/execcasper)
    local_directory="$TMPDIR",  # Use your local directory
    resource_spec="select=1:ncpus=1:mem=23GB",  # Specify resources
    project="ncgd0011",  # Input your project ID here
    walltime="02:00:00",  # Amount of wall time
    interface="ib0",  # Interface to use
)
cluster.scale(6)

/glade/u/home/dcherian/miniconda3/envs/dcpy/lib/python3.8/site-packages/distributed/node.py:160: UserWarning: Port 8787 is already in use.
Perhaps you already have a cluster running?
Hosting the HTTP server on port 46575 instead
  warnings.warn(

client = distributed.Client(cluster)
client

Client

Client-ab60f98f-2082-11ec-920c-3cecef1acbfa

Connection method: Cluster object Cluster type: PBSCluster
Dashboard: /proxy/46575/status

Cluster Info

Stations#

gcmdir = "/glade/campaign/cgd/oce/people/bachman/TPOS_1_20_20_year/OUTPUT/"  # MITgcm output directory
stationdirname = gcmdir

metrics = pump.model.read_metrics(stationdirname)
stations = pump.model.read_stations_20(stationdirname)
enso = pump.obs.make_enso_mask()
stations["enso"] = enso.reindex(time=stations.time.data, method="nearest")

%aimport pump.model

pump.model.rename_metrics(metrics)

<xarray.Dataset>
Dimensions:  (RC: 185, YC: 480, XC: 1500, XG: 1500, YG: 480, RF: 186)
Coordinates:
  * RF       (RF) float64 dask.array<chunksize=(186,), meta=np.ndarray>
Dimensions without coordinates: RC, YC, XC, XG, YG
Data variables: (12/13)
    hFacC    (RC, YC, XC) float64 dask.array<chunksize=(185, 480, 1500), meta=np.ndarray>
    RAC      (YC, XC) float64 dask.array<chunksize=(480, 1500), meta=np.ndarray>
    DXC      (YC, XG) float64 dask.array<chunksize=(480, 1500), meta=np.ndarray>
    DYC      (YG, XC) float64 dask.array<chunksize=(480, 1500), meta=np.ndarray>
    cellvol  (YC, XC, RC) float64 dask.array<chunksize=(480, 1500, 185), meta=np.ndarray>
    rAw      (YC, XG) float32 dask.array<chunksize=(480, 1500), meta=np.ndarray>
    ...       ...
    hFacS    (RC, YG, XC) float32 dask.array<chunksize=(185, 480, 1500), meta=np.ndarray>
    drF      (RC) float32 dask.array<chunksize=(185,), meta=np.ndarray>
    drC      (RF) float32 dask.array<chunksize=(186,), meta=np.ndarray>
    rAs      (YG, XC) float32 dask.array<chunksize=(480, 1500), meta=np.ndarray>
    DXG      (YG, XC) float32 dask.array<chunksize=(480, 1500), meta=np.ndarray>
    DYG      (YC, XG) float32 dask.array<chunksize=(480, 1500), meta=np.ndarray>
Attributes:
    long_name:  cell_height
    units:      m
metrics = metrics.reindex_like(moor, method="nearest").load()

metrics

<xarray.Dataset>
Dimensions:    (depth: 185, latitude: 37, longitude: 4, depth_left: 186, RF: 186)
Coordinates:
  * depth      (depth) float64 -1.25 -3.75 -6.25 ... -5.658e+03 -5.758e+03
  * latitude   (latitude) float64 -3.025 -2.775 -2.525 ... 5.475 5.725 5.975
  * longitude  (longitude) float64 -155.0 -140.0 -125.0 -110.0
    RF         (depth_left) float64 0.0 -2.5 -5.0 ... -5.708e+03 -5.808e+03
Dimensions without coordinates: depth_left
Data variables: (12/14)
    dRF        (depth) float64 -2.5 -2.5 -2.5 -2.5 ... -100.0 -100.0 -100.0
    hFacC      (depth, latitude, longitude) float64 1.0 1.0 1.0 ... 0.0 0.0 0.0
    RAC        (latitude, longitude) float64 3.086e+07 3.086e+07 ... 3.073e+07
    DXC        (latitude, longitude) float64 5.551e+03 5.551e+03 ... 5.529e+03
    DYC        (latitude, longitude) float64 5.559e+03 5.559e+03 ... 5.559e+03
    cellvol    (latitude, longitude, depth) float64 7.715e+07 7.715e+07 ... nan
    ...         ...
    hFacS      (depth, latitude, longitude) float32 1.0 1.0 1.0 ... 0.0 0.0 0.0
    drF        (depth) float32 2.5 2.5 2.5 2.5 2.5 ... 100.0 100.0 100.0 100.0
    drC        (RF) float32 1.25 2.5 2.5 2.5 2.5 ... 100.0 100.0 100.0 50.0
    rAs        (latitude, longitude) float32 3.086e+07 3.086e+07 ... 3.073e+07
    DXG        (latitude, longitude) float32 5.551e+03 5.551e+03 ... 5.529e+03
    DYG        (latitude, longitude) float32 5.559e+03 5.559e+03 ... 5.559e+03
Attributes:
    long_name:  cell_height
    units:      m
moor = stations.isel(longitude=slice(1, -1, 3), latitude=slice(1, -1, 3)).sel(
    time=slice("2007-01-01", "2009-12-31")
)
moor

<xarray.Dataset>
Dimensions:        (depth: 185, time: 26304, longitude: 4, latitude: 37)
Coordinates:
  * depth          (depth) float32 -1.25 -3.75 -6.25 ... -5.658e+03 -5.758e+03
  * latitude       (latitude) float64 -3.025 -2.775 -2.525 ... 5.475 5.725 5.975
  * longitude      (longitude) float64 -155.0 -140.0 -125.0 -110.0
  * time           (time) datetime64[ns] 2007-01-01 ... 2009-12-31T23:00:00
Data variables: (12/21)
    DFrI_TH        (depth, time, longitude, latitude) float32 dask.array<chunksize=(185, 1866, 1, 1), meta=np.ndarray>
    KPPdiffKzT     (depth, time, longitude, latitude) float32 dask.array<chunksize=(185, 1866, 1, 1), meta=np.ndarray>
    KPPg_TH        (depth, time, longitude, latitude) float32 dask.array<chunksize=(185, 1866, 1, 1), meta=np.ndarray>
    KPPhbl         (time, longitude, latitude) float32 dask.array<chunksize=(1866, 1, 1), meta=np.ndarray>
    KPPviscAz      (depth, time, longitude, latitude) float32 dask.array<chunksize=(185, 1866, 1, 1), meta=np.ndarray>
    SSH            (time, longitude, latitude) float32 dask.array<chunksize=(1866, 1, 1), meta=np.ndarray>
    ...             ...
    dens           (depth, time, longitude, latitude) float32 dask.array<chunksize=(185, 1866, 1, 1), meta=np.ndarray>
    mld            (time, longitude, latitude) float32 dask.array<chunksize=(1866, 1, 1), meta=np.ndarray>
    Jq             (depth, time, longitude, latitude) float64 dask.array<chunksize=(185, 1866, 1, 1), meta=np.ndarray>
    dJdz           (depth, time, longitude, latitude) float64 dask.array<chunksize=(185, 1866, 1, 1), meta=np.ndarray>
    dTdt           (depth, time, longitude, latitude) float64 dask.array<chunksize=(185, 1866, 1, 1), meta=np.ndarray>
    enso           (time) <U8 'Neutral' 'Neutral' ... 'El-Nino' 'El-Nino'
Attributes:
    easting:   longitude
    northing:  latitude
    title:     Station profile, index (i,j)=(1201,240)
dim = {"x": "longitude", "y": "latitude", "z": "depth"}
metric = {"x": "DXC", "y": "DYC", "z": "drF"}

der = xr.Dataset()
for axis in ["x", "y", "z"]:
    for var in ["u", "v", "w"]:
        der[f"{var}{axis}"] = stations[var].diff(dim[axis]) / metrics[metric[axis]]
        der[f"{var}{axis}"].attrs["long_name"] = f"${var}_{axis}$"
der = der.reindex_like(moor)
der

/glade/u/home/dcherian/miniconda3/envs/dcpy/lib/python3.8/site-packages/xarray/core/indexing.py:1232: PerformanceWarning: Slicing is producing a large chunk. To accept the large
chunk and silence this warning, set the option
    >>> with dask.config.set(**{'array.slicing.split_large_chunks': False}):
    ...     array[indexer]

To avoid creating the large chunks, set the option
    >>> with dask.config.set(**{'array.slicing.split_large_chunks': True}):
    ...     array[indexer]
  value = value[(slice(None),) * axis + (subkey,)]
/glade/u/home/dcherian/miniconda3/envs/dcpy/lib/python3.8/site-packages/dask/array/core.py:4360: PerformanceWarning: Increasing number of chunks by factor of 37
  result = blockwise(
/glade/u/home/dcherian/miniconda3/envs/dcpy/lib/python3.8/site-packages/xarray/core/indexing.py:1232: PerformanceWarning: Slicing is producing a large chunk. To accept the large
chunk and silence this warning, set the option
    >>> with dask.config.set(**{'array.slicing.split_large_chunks': False}):
    ...     array[indexer]

To avoid creating the large chunks, set the option
    >>> with dask.config.set(**{'array.slicing.split_large_chunks': True}):
    ...     array[indexer]
  value = value[(slice(None),) * axis + (subkey,)]
/glade/u/home/dcherian/miniconda3/envs/dcpy/lib/python3.8/site-packages/dask/array/core.py:4360: PerformanceWarning: Increasing number of chunks by factor of 37
  result = blockwise(
/glade/u/home/dcherian/miniconda3/envs/dcpy/lib/python3.8/site-packages/xarray/core/indexing.py:1232: PerformanceWarning: Slicing is producing a large chunk. To accept the large
chunk and silence this warning, set the option
    >>> with dask.config.set(**{'array.slicing.split_large_chunks': False}):
    ...     array[indexer]

To avoid creating the large chunks, set the option
    >>> with dask.config.set(**{'array.slicing.split_large_chunks': True}):
    ...     array[indexer]
  value = value[(slice(None),) * axis + (subkey,)]
/glade/u/home/dcherian/miniconda3/envs/dcpy/lib/python3.8/site-packages/dask/array/core.py:4360: PerformanceWarning: Increasing number of chunks by factor of 37
  result = blockwise(
/glade/u/home/dcherian/miniconda3/envs/dcpy/lib/python3.8/site-packages/xarray/core/indexing.py:1226: PerformanceWarning: Slicing is producing a large chunk. To accept the large
chunk and silence this warning, set the option
    >>> with dask.config.set(**{'array.slicing.split_large_chunks': False}):
    ...     array[indexer]

To avoid creating the large chunks, set the option
    >>> with dask.config.set(**{'array.slicing.split_large_chunks': True}):
    ...     array[indexer]
  return self.array[key]

<xarray.Dataset>
Dimensions:    (depth: 185, latitude: 37, longitude: 4, time: 26304)
Coordinates:
  * depth      (depth) float64 -1.25 -3.75 -6.25 ... -5.658e+03 -5.758e+03
  * latitude   (latitude) float64 -3.025 -2.775 -2.525 ... 5.475 5.725 5.975
  * longitude  (longitude) float64 -155.0 -140.0 -125.0 -110.0
  * time       (time) datetime64[ns] 2007-01-01 ... 2009-12-31T23:00:00
Data variables:
    ux         (depth, time, longitude, latitude) float64 dask.array<chunksize=(185, 1866, 1, 1), meta=np.ndarray>
    vx         (depth, time, longitude, latitude) float64 dask.array<chunksize=(185, 1866, 1, 1), meta=np.ndarray>
    wx         (depth, time, longitude, latitude) float64 dask.array<chunksize=(185, 1866, 1, 1), meta=np.ndarray>
    uy         (depth, time, longitude, latitude) float64 dask.array<chunksize=(185, 1866, 1, 1), meta=np.ndarray>
    vy         (depth, time, longitude, latitude) float64 dask.array<chunksize=(185, 1866, 1, 1), meta=np.ndarray>
    wy         (depth, time, longitude, latitude) float64 dask.array<chunksize=(185, 1866, 1, 1), meta=np.ndarray>
    uz         (depth, time, longitude, latitude) float32 dask.array<chunksize=(185, 1866, 1, 1), meta=np.ndarray>
    vz         (depth, time, longitude, latitude) float32 dask.array<chunksize=(185, 1866, 1, 1), meta=np.ndarray>
    wz         (depth, time, longitude, latitude) float32 dask.array<chunksize=(185, 1866, 1, 1), meta=np.ndarray>
shred2 = xr.Dataset()
shred2["strvy"] = der.uz**2 * der.vy
shred2["strux"] = der.vz**2 * der.ux
shred2["tiltuy"] = -1 * der.uy * der.vz * der.uz
shred2["tiltvx"] = -1 * der.vx * der.vz * der.uz

shred2["strvy"].attrs["long_name"] = "$u_z^2 * v_y$"
shred2["strux"].attrs["long_name"] = "$v_z^2 * u_x$"
shred2["tiltuy"].attrs["long_name"] = "$-u_y u_z v_z$"
shred2["tiltvx"].attrs["long_name"] = "$- v_x u_z v_z$"

moor = moor.update(shred2).update(der)

subset = (
    moor.sel(latitude=[-2, 0, 2, 4], longitude=-140, method="nearest")
    .sel(depth=slice(-80))
    .mean("depth")
)
subset

/glade/u/home/dcherian/miniconda3/envs/dcpy/lib/python3.8/site-packages/numpy/core/_methods.py:179: RuntimeWarning: invalid value encountered in reduce
  ret = umr_sum(arr, axis, dtype, out, keepdims, where=where)

<xarray.Dataset>
Dimensions:        (time: 26304, latitude: 4)
Coordinates:
  * latitude       (latitude) float64 -2.025 -0.025 1.975 3.975
    longitude      float64 -140.0
  * time           (time) datetime64[ns] 2007-01-01 ... 2009-12-31T23:00:00
Data variables: (12/34)
    DFrI_TH        (time, latitude) float32 dask.array<chunksize=(1866, 1), meta=np.ndarray>
    KPPdiffKzT     (time, latitude) float32 dask.array<chunksize=(1866, 1), meta=np.ndarray>
    KPPg_TH        (time, latitude) float32 dask.array<chunksize=(1866, 1), meta=np.ndarray>
    KPPhbl         (time, latitude) float32 dask.array<chunksize=(1866, 1), meta=np.ndarray>
    KPPviscAz      (time, latitude) float32 dask.array<chunksize=(1866, 1), meta=np.ndarray>
    SSH            (time, latitude) float32 dask.array<chunksize=(1866, 1), meta=np.ndarray>
    ...             ...
    uy             (time, latitude) float64 dask.array<chunksize=(1866, 1), meta=np.ndarray>
    vy             (time, latitude) float64 dask.array<chunksize=(1866, 1), meta=np.ndarray>
    wy             (time, latitude) float64 dask.array<chunksize=(1866, 1), meta=np.ndarray>
    uz             (time, latitude) float32 dask.array<chunksize=(1866, 1), meta=np.ndarray>
    vz             (time, latitude) float32 dask.array<chunksize=(1866, 1), meta=np.ndarray>
    wz             (time, latitude) float32 dask.array<chunksize=(1866, 1), meta=np.ndarray>
Attributes:
    easting:   longitude
    northing:  latitude
    title:     Station profile, index (i,j)=(1201,240)
%aimport dask_groupby.core
%aimport dask_groupby.xarray
from dask_groupby.xarray import resample_reduce

%aimport dcpy.ts

daily = resample_reduce(subset.resample(time="D"), func="mean").compute()

dcpy.ts.PlotCoherence(sub.vx, sub.vz, unwrap=True);
_images/6fcc01524125e41dabbf72fe1f5a3674c5e06102a4d9558703c1762440c53859.png

Need indicators for phenomena#

sub = daily.sel(latitude=0, method="nearest")
_, ax = dcpy.ts.PlotSpectrum(sub.v)
dcpy.ts.PlotSpectrum(sub.vx * 1e5, ax=ax)
dcpy.ts.PlotSpectrum(sub.vz * 100, ax=ax)

([<matplotlib.lines.Line2D at 0x2ac4e9310ca0>],
 <AxesSubplot:title={'center':' $v_y$'}, xlabel='Frequency [cpd]', ylabel='PSD'>)
_images/5720dc78bf578b6b7be65c64a460b3349165ac842ad1b05606f57edb48884d8d.png _images/8d32496a966299a1054a5874c82e427a8fe441c60377009f791aec16cedb4bde.png

Divergence terms#

What is happening at 7 day periods? \(v_y \sim w_z\)?

sub = daily.sel(latitude=0, method="nearest")
_, ax = dcpy.ts.PlotSpectrum(sub.ux)
dcpy.ts.PlotSpectrum(sub.vy, ax=ax)
dcpy.ts.PlotSpectrum(sub.wz, ax=ax)

([<matplotlib.lines.Line2D at 0x2ac4e84ef610>],
 <AxesSubplot:title={'center':' $w_z$'}, xlabel='Frequency [cpd]', ylabel='PSD'>)
_images/a01b6a1784f0b6276dc42c60050208645f67ee446a37aa8ad69e98a5e4736477.png

Holmes et al: vortex stretching term \(u_z v_y\) intensifies \(-u_z\)#

sub = daily.sel(latitude=0, method="nearest")
term = sub.strvy
_, ax = dcpy.ts.PlotCoherence(term, sub.v)
dcpy.ts.PlotCoherence(term, sub.vz, ax=ax, ploty0=False, unwrap=True)
dcpy.ts.PlotCoherence(term, sub.vx, ax=ax, ploty0=False, unwrap=True);
_images/b19c2ec169c81ac3b6c702a6e6dbdf8346023da3fbdf7575f8e31fadcd83fe2a.png

Cherian et al: vortex tilting term \(-u_y u_z\) rotates \(u_z\) to create \(v_z\) off the equator#

sub = daily.sel(latitude=3, method="nearest")

term = sub.tiltuy
_, ax = dcpy.ts.PlotCoherence(term, sub.v)
dcpy.ts.PlotCoherence(
    term,
    daily.vz.sel(latitude=0, method="nearest"),
    ax=ax,
    ploty0=False,
    unwrap=True,
)
dcpy.ts.PlotCoherence(
    term,
    daily.vx.sel(latitude=0, method="nearest"),
    ax=ax,
    ploty0=False,
    unwrap=True,
);
_images/0c369ff51f04ee8a3511d4f6a6fa2ba349f24b6d6ec6b8ab3af56a62d419b4e1.png 
sub = daily.sel(latitude=0, method="nearest")

term = sub.vy * sub.uz
_, ax = dcpy.ts.PlotCoherence(term, sub.v)
dcpy.ts.PlotCoherence(
    term,
    daily.vz.sel(latitude=0, method="nearest"),
    ax=ax,
    ploty0=False,
    unwrap=True,
)
dcpy.ts.PlotCoherence(
    term,
    daily.vx.sel(latitude=0, method="nearest"),
    ax=ax,
    ploty0=False,
    unwrap=True,
);
_images/56c9c7a5ffab0ccfe2b819f44a83e9eb68257a262f5c9a5703d2546fe04fecb2.png

Checking all terms#

sub = daily.sel(latitude=0, method="nearest")


def plot_spec(ts, ax):
    dcpy.ts.PlotSpectrum(ts, ax=ax, label=ts.attrs["long_name"])


f, ax = plt.subplots(1, 1)
for term in [sub.strux, sub.strvy, sub.tiltuy, sub.tiltvx]:
    plot_spec(term, ax)
ax.legend()

<matplotlib.legend.Legend at 0x2ac4f0364e20>
_images/fc64bd3b07a0985ec17fe596a7ec466bf37691c2c93d71ebaac2de39005fe1ea.png 
sub = daily.sel(latitude=3.5, method="nearest")


def plot_spec(ts, ax):
    dcpy.ts.PlotSpectrum(ts, ax=ax, label=ts.attrs["long_name"])


f, ax = plt.subplots(1, 1)
for term in [sub.strux, sub.strvy, sub.tiltuy, sub.tiltvx]:
    plot_spec(term, ax)
ax.legend()

<matplotlib.legend.Legend at 0x2ac4e633d310>
_images/f05afccf59493ae3fd5d5e45b76f074ce4a5d676cd0c25e6f1baf285db8b548b.png 
(
    daily[["strux", "strvy", "tiltuy", "tiltvx"]]
    .sel(latitude=4, method="nearest")
    .to_array(dim="variable")
    .plot(hue="variable")
)

[<matplotlib.lines.Line2D at 0x2ac4e9cfce80>,
 <matplotlib.lines.Line2D at 0x2ac4e9cfcf70>,
 <matplotlib.lines.Line2D at 0x2ac4e9cfcf40>,
 <matplotlib.lines.Line2D at 0x2ac4e9d060d0>]
_images/5f424d804ef9ece69dae1f9d1f64837cf797dbf20fe115320fe43b89495be301.png 
station_grid = xgcm.Grid(
    stations.merge(station_metrics),
    coords={
        "X": {"center": "longitude"},
        "Y": {"center": "latitude"},
        "Z": {"center": "depth"},
    },
    metrics={"X": ("DXC",), "Y": ("DYC",), "Z": ("drC",), ("X", "Y"): ("RAC",)},
    periodic=False,
    boundary="fill",
    fill_value=np.nan,
)

Full domain#

Having lots of trouble getting this to work.

ds, metrics, grid = pump.model.read_mitgcm_20_year(
    state=True,
    start="2008-01-01",
    stop="2008-12-31",
    chunks={"latitude": -1, "longitude": -1},
)
ds["b"] = -9.81 / 1035 * ds.dens
ds

3288 3653

<xarray.Dataset>
Dimensions:          (RF: 136, YG: 400, XG: 1420, YC: 400, XC: 1420, time: 366, RC: 136)
Coordinates:
  * RF               (RF) float64 0.0 -2.5 -5.0 -7.5 ... -797.0 -851.7 -911.6
  * YG               (YG) float64 -10.0 -9.95 -9.9 -9.85 ... 9.8 9.85 9.9 9.95
  * XG               (XG) float64 -168.0 -168.0 -167.9 ... -97.18 -97.12 -97.08
  * YC               (YC) float64 -9.975 -9.925 -9.875 ... 9.875 9.925 9.975
  * XC               (XC) float64 -168.0 -167.9 -167.9 ... -97.15 -97.1 -97.05
  * time             (time) datetime64[ns] 2008-01-01 2008-01-02 ... 2008-12-31
  * RC               (RC) float64 -1.25 -3.75 -6.25 ... -824.4 -881.7 -944.4
Data variables: (12/14)
    theta            (time, RC, YC, XC) float32 dask.array<chunksize=(1, 136, 400, 1420), meta=np.ndarray>
    u                (time, RC, YC, XG) float32 dask.array<chunksize=(1, 136, 400, 1420), meta=np.ndarray>
    v                (time, RC, YG, XC) float32 dask.array<chunksize=(1, 136, 400, 1420), meta=np.ndarray>
    w                (time, RF, YC, XC) float32 dask.array<chunksize=(1, 136, 400, 1420), meta=np.ndarray>
    salt             (time, RC, YC, XC) float32 dask.array<chunksize=(1, 136, 400, 1420), meta=np.ndarray>
    KPP_diffusivity  (time, RF, YC, XC) float32 dask.array<chunksize=(1, 136, 400, 1420), meta=np.ndarray>
    ...               ...
    S2               (time, RF, YC, XC) float32 dask.array<chunksize=(1, 1, 399, 1419), meta=np.ndarray>
    Ri               (time, RF, YC, XC) float32 dask.array<chunksize=(1, 1, 399, 1419), meta=np.ndarray>
    mld              (time, YC, XC) float64 dask.array<chunksize=(1, 400, 1420), meta=np.ndarray>
    dcl_base         (time, YC, XC) float64 dask.array<chunksize=(1, 399, 1419), meta=np.ndarray>
    dcl              (time, YC, XC) float64 dask.array<chunksize=(1, 399, 1419), meta=np.ndarray>
    b                (time, RC, YC, XC) float32 dask.array<chunksize=(1, 136, 400, 1420), meta=np.ndarray>
Attributes:
    title:              daily snapshot from 1/20 degree Equatorial Pacific MI...
    easting:            longitude
    northing:           latitude
    field_julian_date:  508416
    julian_day_unit:    days since 1950-01-01 00:00:00
ds["dcl"], ds["dcl_base"] = dask.persist(ds.dcl, ds.dcl_base)

der = xr.Dataset()
for axis in ["x", "y", "z"]:
    for field in ["u", "v", "w", "b"]:
        der[f"{field}{axis}"] = grid.derivative(ds[field], axis=axis.upper())
der

/glade/u/home/dcherian/python/xgcm/xgcm/grid.py:1384: UserWarning: Metric at ('time', 'RC', 'YC', 'XC') being interpolated from metrics at dimensions ('YG', 'XC'). Boundary value set to 'extend'.
  warnings.warn(
/glade/u/home/dcherian/python/xgcm/xgcm/grid.py:1384: UserWarning: Metric at ('time', 'RC', 'YG', 'XG') being interpolated from metrics at dimensions ('YG', 'XC'). Boundary value set to 'extend'.
  warnings.warn(
/glade/u/home/dcherian/python/xgcm/xgcm/grid.py:1384: UserWarning: Metric at ('time', 'RC', 'YG', 'XG') being interpolated from metrics at dimensions ('YG', 'XC'). Boundary value set to 'extend'.
  warnings.warn(
/glade/u/home/dcherian/python/xgcm/xgcm/grid.py:1384: UserWarning: Metric at ('time', 'RC', 'YC', 'XC') being interpolated from metrics at dimensions ('YG', 'XC'). Boundary value set to 'extend'.
  warnings.warn(

<xarray.Dataset>
Dimensions:  (time: 366, RC: 136, YC: 400, XC: 1420, YG: 400, XG: 1420, RF: 136)
Coordinates:
  * time     (time) datetime64[ns] 2008-01-01 2008-01-02 ... 2008-12-31
  * RC       (RC) float64 -1.25 -3.75 -6.25 -8.75 ... -824.4 -881.7 -944.4
  * YC       (YC) float64 -9.975 -9.925 -9.875 -9.825 ... 9.875 9.925 9.975
  * XC       (XC) float64 -168.0 -167.9 -167.9 -167.9 ... -97.15 -97.1 -97.05
  * YG       (YG) float64 -10.0 -9.95 -9.9 -9.85 -9.8 ... 9.75 9.8 9.85 9.9 9.95
  * XG       (XG) float64 -168.0 -168.0 -167.9 -167.9 ... -97.18 -97.12 -97.08
  * RF       (RF) float64 0.0 -2.5 -5.0 -7.5 ... -747.1 -797.0 -851.7 -911.6
Data variables:
    ux       (time, RC, YC, XC) float32 dask.array<chunksize=(1, 136, 400, 1419), meta=np.ndarray>
    vx       (time, RC, YG, XG) float32 dask.array<chunksize=(1, 136, 400, 1), meta=np.ndarray>
    wx       (time, RF, YC, XG) float64 dask.array<chunksize=(1, 136, 400, 1), meta=np.ndarray>
    bx       (time, RC, YC, XG) float64 dask.array<chunksize=(1, 136, 400, 1), meta=np.ndarray>
    uy       (time, RC, YG, XG) float64 dask.array<chunksize=(1, 136, 1, 1420), meta=np.ndarray>
    vy       (time, RC, YC, XC) float64 dask.array<chunksize=(1, 136, 399, 1420), meta=np.ndarray>
    wy       (time, RF, YG, XC) float64 dask.array<chunksize=(1, 136, 1, 1420), meta=np.ndarray>
    by       (time, RC, YG, XC) float64 dask.array<chunksize=(1, 136, 1, 1420), meta=np.ndarray>
    uz       (time, RF, YC, XG) float32 dask.array<chunksize=(1, 1, 400, 1420), meta=np.ndarray>
    vz       (time, RF, YG, XC) float32 dask.array<chunksize=(1, 1, 400, 1420), meta=np.ndarray>
    wz       (time, RC, YC, XC) float32 dask.array<chunksize=(1, 135, 400, 1420), meta=np.ndarray>
    bz       (time, RF, YC, XC) float32 dask.array<chunksize=(1, 1, 400, 1420), meta=np.ndarray>
shred2 = xr.Dataset()
shred2["strvy"] = grid.interp_like(der.uz**2, der.vy) * der.vy
shred2["strux"] = grid.interp_like(der.vz**2, der.ux) * der.ux
shred2["str"] = shred2.strvy + shred2.strux

shred2["tilt"] = (
    -grid.interp_like(der.uy + der.vx, shred2.str)
    * grid.interp_like(der.uz, shred2.str)
    * grid.interp_like(der.vz, shred2.str)
)
shred2

<xarray.Dataset>
Dimensions:  (time: 366, RC: 136, YC: 400, XC: 1420)
Coordinates:
  * time     (time) datetime64[ns] 2008-01-01 2008-01-02 ... 2008-12-31
  * RC       (RC) float64 -1.25 -3.75 -6.25 -8.75 ... -824.4 -881.7 -944.4
  * YC       (YC) float64 -9.975 -9.925 -9.875 -9.825 ... 9.875 9.925 9.975
  * XC       (XC) float64 -168.0 -167.9 -167.9 -167.9 ... -97.15 -97.1 -97.05
Data variables:
    strvy    (time, RC, YC, XC) float64 dask.array<chunksize=(1, 1, 399, 1419), meta=np.ndarray>
    strux    (time, RC, YC, XC) float32 dask.array<chunksize=(1, 1, 399, 1419), meta=np.ndarray>
    str      (time, RC, YC, XC) float64 dask.array<chunksize=(1, 1, 399, 1419), meta=np.ndarray>
    tilt     (time, RC, YC, XC) float64 dask.array<chunksize=(1, 1, 1, 1), meta=np.ndarray>
dclsub = ds.dcl
dclsub = dclsub.where(dclsub > 10)
dclsub.isel(time=0).plot(robust=True)

<matplotlib.collections.QuadMesh at 0x2b3464ff8a00>
_images/0af2e78ac08e6cc1b8c34ebbb3042a75fe4e2d2d3d76cb73c3282501dbdad8b8.png 
shred2140 = dcpy.dask.map_copy(
    shred2.strvy.sel(YC=0, XC=-140, method="nearest")
).persist()

dclmask = (ds.dcl > 10) & (ds.RC > ds.dcl_base)
dclmask

<xarray.DataArray (time: 366, YC: 400, XC: 1420, RC: 136)>
dask.array<and_, shape=(366, 400, 1420, 136), dtype=bool, chunksize=(1, 399, 1419, 136), chunktype=numpy.ndarray>
Coordinates:
  * YC       (YC) float64 -9.975 -9.925 -9.875 -9.825 ... 9.875 9.925 9.975
  * XC       (XC) float64 -168.0 -167.9 -167.9 -167.9 ... -97.15 -97.1 -97.05
  * time     (time) datetime64[ns] 2008-01-01 2008-01-02 ... 2008-12-31
  * RC       (RC) float64 -1.25 -3.75 -6.25 -8.75 ... -824.4 -881.7 -944.4
Attributes:
    long_name:    DCL
    units:        m
    description:  Interpolate density to 1m grid. Search for max depth where ...
mivol = grid.integrate(dclmask, axis=("X", "Y", "Z"))
mivol

/glade/u/home/dcherian/python/xgcm/xgcm/grid.py:1408: UserWarning: Metric at ('time', 'YC', 'XC', 'RC') being interpolated from metrics at dimensions [('YC', 'XG'), ('RC',)]. Boundary value set to 'extend'.
  warnings.warn(
/glade/u/home/dcherian/python/xgcm/xgcm/grid.py:1408: UserWarning: Metric at ('time', 'YC', 'XC', 'RC') being interpolated from metrics at dimensions [('YC', 'XG'), ('RF',)]. Boundary value set to 'extend'.
  warnings.warn(
/glade/u/home/dcherian/python/xgcm/xgcm/grid.py:1408: UserWarning: Metric at ('time', 'YC', 'XC', 'RC') being interpolated from metrics at dimensions [('YG', 'XC'), ('RC',)]. Boundary value set to 'extend'.
  warnings.warn(
/glade/u/home/dcherian/python/xgcm/xgcm/grid.py:1408: UserWarning: Metric at ('time', 'YC', 'XC', 'RC') being interpolated from metrics at dimensions [('YG', 'XC'), ('RF',)]. Boundary value set to 'extend'.
  warnings.warn(

<xarray.DataArray (time: 366)>
dask.array<sum-aggregate, shape=(366,), dtype=float32, chunksize=(1,), chunktype=numpy.ndarray>
Coordinates:
  * time     (time) datetime64[ns] 2008-01-01 2008-01-02 ... 2008-12-31
Attributes:
    long_name:    DCL
    units:        m
    description:  Interpolate density to 1m grid. Search for max depth where ...
intshred2 = grid.integrate(shred2.where(dclmask), axis=("X", "Y", "Z"))
intshred2

/glade/u/home/dcherian/python/xgcm/xgcm/grid.py:1408: UserWarning: Metric at Frozen({'time': 366, 'RC': 136, 'YC': 400, 'XC': 1420}) being interpolated from metrics at dimensions [('YC', 'XG'), ('RC',)]. Boundary value set to 'extend'.
  warnings.warn(
/glade/u/home/dcherian/python/xgcm/xgcm/grid.py:1408: UserWarning: Metric at Frozen({'time': 366, 'RC': 136, 'YC': 400, 'XC': 1420}) being interpolated from metrics at dimensions [('YC', 'XG'), ('RF',)]. Boundary value set to 'extend'.
  warnings.warn(
/glade/u/home/dcherian/python/xgcm/xgcm/grid.py:1408: UserWarning: Metric at Frozen({'time': 366, 'RC': 136, 'YC': 400, 'XC': 1420}) being interpolated from metrics at dimensions [('YG', 'XC'), ('RC',)]. Boundary value set to 'extend'.
  warnings.warn(
/glade/u/home/dcherian/python/xgcm/xgcm/grid.py:1408: UserWarning: Metric at Frozen({'time': 366, 'RC': 136, 'YC': 400, 'XC': 1420}) being interpolated from metrics at dimensions [('YG', 'XC'), ('RF',)]. Boundary value set to 'extend'.
  warnings.warn(

<xarray.Dataset>
Dimensions:  (time: 366)
Coordinates:
  * time     (time) datetime64[ns] 2008-01-01 2008-01-02 ... 2008-12-31
Data variables:
    strvy    (time) float64 dask.array<chunksize=(1,), meta=np.ndarray>
    strux    (time) float32 dask.array<chunksize=(1,), meta=np.ndarray>
    str      (time) float64 dask.array<chunksize=(1,), meta=np.ndarray>
    tilt     (time) float64 dask.array<chunksize=(1,), meta=np.ndarray>
mivol, intshred2 = dask.persist(mivol, intshred2)

cluster.scale(24)