TIW vertical velocity in MITgcm simulations#

TODO:

  1. Conditional average by TIW KE?

  2. look more closely at downwelling in 2008 OND

  3. regrid to ML-relative vertical coordinate; not necessary MLD seems quite small!

  4. average 2005,2006; and then 2008, 2009 because they look similar

  5. better season definition: based on TIWKE?

Hide code cell content 
# %load_ext %autoreload
%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))


print(dask.__version__)
print(distributed.__version__)
print(xr.__version__)


xr.DataArray([1.0])

2021.07.2
2021.07.2
0.19.0

<xarray.DataArray (dim_0: 1)>
array([1.])
Dimensions without coordinates: dim_0
Hide code cell content 
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)
Hide code cell content 
cluster.scale(24)
Hide code cell content 
client = distributed.Client(cluster)
client

Client

Client-1880c442-116b-11ec-b9b4-3cecef1b11de

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

Cluster Info

Means#

ds, metrics = pump.model.read_mitgcm_20_year(
    state=True, start="2008-01-01", stop="2017-12-31"
)
ds

seasonalmean = ds.resample(time="QS").mean()
seasonalmean

<xarray.Dataset>
Dimensions:          (time: 68, RF: 136, YG: 400, XG: 1420, YC: 400, XC: 1420, RC: 136)
Coordinates:
  * time             (time) datetime64[ns] 2001-01-01 2001-04-01 ... 2017-10-01
  * 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
  * RC               (RC) float64 -1.25 -3.75 -6.25 ... -824.4 -881.7 -944.4
Data variables:
    theta            (time, RC, YC, XC) float32 dask.array<chunksize=(1, 136, 400, 500), meta=np.ndarray>
    u                (time, RC, YC, XG) float32 dask.array<chunksize=(1, 136, 400, 500), meta=np.ndarray>
    v                (time, RC, YG, XC) float32 dask.array<chunksize=(1, 136, 400, 500), meta=np.ndarray>
    w                (time, RF, YC, XC) float32 dask.array<chunksize=(1, 136, 400, 500), meta=np.ndarray>
    salt             (time, RC, YC, XC) float32 dask.array<chunksize=(1, 136, 400, 500), meta=np.ndarray>
    KPP_diffusivity  (time, RF, YC, XC) float32 dask.array<chunksize=(1, 136, 400, 500), meta=np.ndarray>
    dens             (time, RC, YC, XC) float32 dask.array<chunksize=(1, 136, 400, 500), meta=np.ndarray>
Attributes:
    title:              daily snapshot from 1/20 degree Equatorial Pacific MI...
    easting:            longitude
    northing:           latitude
    field_julian_date:  447072
    julian_day_unit:    days since 1950-01-01 00:00:00
dcpy.dask.batch_to_zarr(
    seasonalmean[["u", "theta"]],
    file="/glade/work/dcherian/pump/zarrs/seasonalutheta.zarr",
    batch_size=12,
    dim="time",
)

100%|██████████| 6/6 [20:23<00:00, 203.96s/it]

cluster.scale(12)

v = dcpy.dask.map_copy(ds.v.sel(RC=slice(-80)).sel(XC=-140, YG=0, method="nearest"))
v = v.load()
v.RC.attrs["axis"] = "Z"

v.to_dataset().to_zarr("/glade/work/dcherian/pump/zarrs/v-0-140-80m.zarr")

<xarray.backends.zarr.ZarrStore at 0x2ba17e0f3ba0>

TIWKE#

make TIW KE proxy for later conditional averaging

def tiw_avg_filter_v(v):
    import xfilter

    v = xfilter.bandpass(
        v.cf.sel(Z=slice(-10, -80)).cf.mean("Z"),
        coord="time",
        freq=[1 / 10.0, 1 / 50.0],
        cycles_per="D",
        method="pad",
        num_discard=0,
        order=3,
    )

    if v.count() == 0:
        raise ValueError("No good data in filtered depth-averaged v.")

    v.attrs["long_name"] = "v: (10, 80m) avg, 10d-40d bandpass"

    return v


# tiwv = tiw_avg_filter_v(v)

dcpy.ts.PlotSpectrum(v)

([<matplotlib.lines.Line2D at 0x2ba1f7dc9460>,
  <matplotlib.lines.Line2D at 0x2ba1f2f34b50>,
  <matplotlib.lines.Line2D at 0x2ba1f6d644f0>,
  <matplotlib.lines.Line2D at 0x2ba1f2f34e50>,
  <matplotlib.lines.Line2D at 0x2ba1f2f34310>,
  <matplotlib.lines.Line2D at 0x2ba1fc7fc940>,
  <matplotlib.lines.Line2D at 0x2ba1fd14f250>,
  <matplotlib.lines.Line2D at 0x2ba1fcea8610>,
  <matplotlib.lines.Line2D at 0x2ba20e8b7040>,
  <matplotlib.lines.Line2D at 0x2ba1fc26d940>,
  <matplotlib.lines.Line2D at 0x2ba181a6c190>,
  <matplotlib.lines.Line2D at 0x2ba1fe1d8ee0>,
  <matplotlib.lines.Line2D at 0x2ba0d09da340>,
  <matplotlib.lines.Line2D at 0x2ba2122b1a00>,
  <matplotlib.lines.Line2D at 0x2ba210695460>,
  <matplotlib.lines.Line2D at 0x2ba20f1e6190>,
  <matplotlib.lines.Line2D at 0x2ba1f69d09d0>,
  <matplotlib.lines.Line2D at 0x2ba1fca98460>,
  <matplotlib.lines.Line2D at 0x2ba1fe122100>,
  <matplotlib.lines.Line2D at 0x2ba0d09da2e0>,
  <matplotlib.lines.Line2D at 0x2ba1f7c02850>,
  <matplotlib.lines.Line2D at 0x2ba2122b1820>,
  <matplotlib.lines.Line2D at 0x2ba1fcea87f0>,
  <matplotlib.lines.Line2D at 0x2ba1fcb76700>,
  <matplotlib.lines.Line2D at 0x2ba1f68f3a30>,
  <matplotlib.lines.Line2D at 0x2ba1fc7d4a00>,
  <matplotlib.lines.Line2D at 0x2ba1a7db2310>,
  <matplotlib.lines.Line2D at 0x2ba20c1c7c70>,
  <matplotlib.lines.Line2D at 0x2ba1fc7235b0>,
  <matplotlib.lines.Line2D at 0x2ba1fc26d640>,
  <matplotlib.lines.Line2D at 0x2ba1fc7238b0>,
  <matplotlib.lines.Line2D at 0x2ba1fc7d4370>],
 <AxesSubplot:title={'center':' v'}, xlabel='Wavelength/2π', ylabel='PSD'>)
_images/9a9fcffb56f75abc32336259a2ffae2550243765665e25cdcddfff183074d359.png 
tiwke = dask_groupby.xarray.resample_reduce(
    (tiwv**2 / 2).resample(time="M"), func="mean"
)
tiwke.plot()

[<matplotlib.lines.Line2D at 0x2ba2138bb370>]
_images/386b80466af02f6e468535fec55d1334823016a44755d76a8f187ea9d07dddde.png

Annual means#

wmean = xr.open_zarr("meanw.zarr")

wmean.sel(XC=-140, method="nearest").sel(YC=slice(-5, 5), RF=slice(-200)).w.plot(
    col="time",
    col_wrap=4,
    y="RF",
    x="YC",
    robust=True,
    cmap="RdBu_r",
    cbar_kwargs={"orientation": "horizontal"},
)

<xarray.plot.facetgrid.FacetGrid at 0x2ba0d082a5e0>
_images/2f2d1bbf8fb565248a0dd6ff0ce9b3fe25aa79c8ee3ed54cea70b6f6b2986869.png 
wmean.sel(XC=-140, method="nearest").sel(YC=slice(-5, 5), RF=slice(-200)).w.plot(
    col="time",
    col_wrap=3,
    y="RF",
    x="YC",
    robust=True,
    cmap="RdBu_r",
    cbar_kwargs={"orientation": "horizontal"},
)

<xarray.plot.facetgrid.FacetGrid at 0x2b8a6f5464f0>
_images/114236fd9f5b9e8e17b194b5c37cd4f60bff961b834d5d56113c77c6547e07e3.png

Seasonal mean#

Hide code cell content 
seasonalmean = xr.open_mfdataset("/glade/work/dcherian/pump/zarrs/seasonal*.zarr")
seasonalmean

/glade/u/home/dcherian/miniconda3/envs/dcpy/lib/python3.8/site-packages/xarray/backends/plugins.py:110: RuntimeWarning: 'netcdf4' fails while guessing
  warnings.warn(f"{engine!r} fails while guessing", RuntimeWarning)
/glade/u/home/dcherian/miniconda3/envs/dcpy/lib/python3.8/site-packages/xarray/backends/plugins.py:110: RuntimeWarning: 'h5netcdf' fails while guessing
  warnings.warn(f"{engine!r} fails while guessing", RuntimeWarning)
/glade/u/home/dcherian/miniconda3/envs/dcpy/lib/python3.8/site-packages/xarray/backends/plugins.py:110: RuntimeWarning: 'scipy' fails while guessing
  warnings.warn(f"{engine!r} fails while guessing", RuntimeWarning)
/glade/u/home/dcherian/miniconda3/envs/dcpy/lib/python3.8/site-packages/xarray/backends/plugins.py:110: RuntimeWarning: 'netcdf4' fails while guessing
  warnings.warn(f"{engine!r} fails while guessing", RuntimeWarning)
/glade/u/home/dcherian/miniconda3/envs/dcpy/lib/python3.8/site-packages/xarray/backends/plugins.py:110: RuntimeWarning: 'h5netcdf' fails while guessing
  warnings.warn(f"{engine!r} fails while guessing", RuntimeWarning)
/glade/u/home/dcherian/miniconda3/envs/dcpy/lib/python3.8/site-packages/xarray/backends/plugins.py:110: RuntimeWarning: 'scipy' fails while guessing
  warnings.warn(f"{engine!r} fails while guessing", RuntimeWarning)

<xarray.Dataset>
Dimensions:  (RC: 136, XC: 1420, XG: 1420, YC: 400, time: 68, RF: 136)
Coordinates:
  * RC       (RC) float64 -1.25 -3.75 -6.25 -8.75 ... -824.4 -881.7 -944.4
  * XC       (XC) float64 -168.0 -167.9 -167.9 -167.9 ... -97.15 -97.1 -97.05
  * XG       (XG) float64 -168.0 -168.0 -167.9 -167.9 ... -97.18 -97.12 -97.08
  * YC       (YC) float64 -9.975 -9.925 -9.875 -9.825 ... 9.875 9.925 9.975
  * time     (time) datetime64[ns] 2001-01-01 2001-04-01 ... 2017-10-01
  * RF       (RF) float64 0.0 -2.5 -5.0 -7.5 ... -747.1 -797.0 -851.7 -911.6
Data variables:
    theta    (time, RC, YC, XC) float32 dask.array<chunksize=(1, 136, 400, 500), meta=np.ndarray>
    u        (time, RC, YC, XG) float32 dask.array<chunksize=(1, 136, 400, 500), meta=np.ndarray>
    w        (time, RF, YC, XC) float32 dask.array<chunksize=(1, 136, 400, 500), meta=np.ndarray>
Attributes:
    easting:            longitude
    field_julian_date:  447072
    julian_day_unit:    days since 1950-01-01 00:00:00
    northing:           latitude
    title:              daily snapshot from 1/20 degree Equatorial Pacific MI...
seasonal140 = (
    seasonalmean.groupby("time.month")
    .mean()
    .sel(YC=slice(-5, 5), RF=slice(-250))
    .sel(XC=-140, XG=-140, method="nearest")
).load()

fg = seasonal140.w.plot(
    col="month",
    col_wrap=4,
    y="RF",
    x="YC",
    robust=True,
    cmap="RdBu_r",
    cbar_kwargs={"orientation": "horizontal"},
)

fg.data = seasonal140.u
fg.map_dataarray(
    xr.plot.contour, levels=11, x="YC", y="RC", colors="k", add_colorbar=False
)
fg.data = seasonal140.theta
fg.map_dataarray(
    xr.plot.contour, levels=11, x="YC", y="RC", colors="g", add_colorbar=False
)

<xarray.plot.facetgrid.FacetGrid at 0x2ba1723c5820>
_images/7b7df6fe276165d7da815e3bd1fb3ebe4fbe0656cab7508118ac62f0a7e8b154.png

Seasonal mean by year#

seasonalmean.sel(time=slice("2005", "2009")).sel(XC=-140, method="nearest").sel(
    YC=slice(-5, 5), RF=slice(-250)
).w.plot(
    col="time",
    col_wrap=4,
    y="RF",
    x="YC",
    robust=True,
    cmap="RdBu_r",
    cbar_kwargs={"orientation": "horizontal"},
)

<xarray.plot.facetgrid.FacetGrid at 0x2ba17b6cb2b0>
_images/95589167ff3cb7e85e2844cf2947af34c0be078fc7390b00b49ee595c2048b01.png

Seasonal mean by ENSO phase#

TODO I am picking the enso phase in the first month of the season. I could do this in a better way…

Hide code cell content 
enso_phase = pump.obs.make_enso_mask()
seasonalmean.coords["enso_phase"] = (
    enso_phase.resample(time="QS").first().reindex(time=seasonalmean.time)
)
Hide code cell content 
import dask_groupby.xarray

seasonal_phase_mean = dask_groupby.xarray.xarray_reduce(
    seasonalmean.sel(XC=[-140], XG=[-140], method="nearest"),
    seasonalmean.time.dt.month,
    seasonalmean.enso_phase,
    func="mean",
    fill_value=np.nan,
)
Hide code cell content 
seasonal_phase_mean.load()

<xarray.Dataset>
Dimensions:     (YC: 400, RC: 136, XC: 1, month: 4, enso_phase: 3, XG: 1, RF: 136)
Coordinates:
  * RC          (RC) float64 -1.25 -3.75 -6.25 -8.75 ... -824.4 -881.7 -944.4
  * XC          (XC) float64 -140.0
  * XG          (XG) float64 -140.0
  * YC          (YC) float64 -9.975 -9.925 -9.875 -9.825 ... 9.875 9.925 9.975
  * RF          (RF) float64 0.0 -2.5 -5.0 -7.5 ... -747.1 -797.0 -851.7 -911.6
  * month       (month) int64 1 4 7 10
  * enso_phase  (enso_phase) object 'La-Nina' 'Neutral' 'El-Nino'
Data variables:
    theta       (YC, RC, XC, month, enso_phase) float64 27.42 27.86 ... 4.675
    u           (YC, RC, XG, month, enso_phase) float64 -0.09301 ... -0.01501
    w           (YC, RF, XC, month, enso_phase) float64 -5.711e-09 ... 7.112e-06
Attributes:
    easting:            longitude
    field_julian_date:  447072
    julian_day_unit:    days since 1950-01-01 00:00:00
    northing:           latitude
    title:              daily snapshot from 1/20 degree Equatorial Pacific MI...

  1. There is definitely a northward “bias” in upwelling

  2. The split is more obvious in La-Nina JFM, OND seasons?

Hide code cell source 
fg = (
    seasonal_phase_mean.sel(YC=slice(-7, 7), RF=slice(-250))
    .sel(XC=-140, method="nearest")
    .w.plot(
        col="month",
        row="enso_phase",
        y="RF",
        x="YC",
        robust=True,
        cmap="RdBu_r",
        cbar_kwargs={"orientation": "horizontal"},
    )
)

fg.data = (
    seasonal_phase_mean.sel(YC=slice(-7, 7), RF=slice(-250))
    .sel(XG=-140, method="nearest")
    .u
)

fg.map_dataarray(
    xr.plot.contour,
    y="RC",
    x="YC",
    levels=21,
    colors="k",
    linewidths=0.5,
    add_colorbar=False,
)

<xarray.plot.facetgrid.FacetGrid at 0x2ba17b4b29a0>
_images/aa61be9311d5824c229301a34e2534f2bd3b39c316c7aaade4840f71a783d66f.png

Lets check the OND season for El-Nino years

seasonalmean.sel(XC=-140, method="nearest").sel(
    time=seasonalmean.time.dt.month == 10
).query(time="enso_phase == 'El-Nino'").w.sel(RF=slice(-200), YC=slice(-5, 5)).plot(
    col="time",
    y="RF",
    x="YC",
    robust=True,
    cmap="RdBu_r",
    cbar_kwargs={"orientation": "vertical"},
)

<xarray.plot.facetgrid.FacetGrid at 0x2ba182b72fd0>
_images/e188c047322e15ef36cec1828da52017e0665ec338ef5aa9a0780c2973469fed.png

2005 vs 2008 OND at 140W#

2008 has a more pronounced downwelling signal at the equator unlike 2005. But turns out they both have a bunch of TIWs

TODO extend this for 2005, 2006, 2007, 2008.

read data#

Hide code cell content 
ond05, metrics, grid = pump.model.read_mitgcm_20_year(
    state=True,
    start="2005-10-01",
    stop="2005-12-31",  # chunks={"RC": 10, "RF": 10},
)
ond05

2466 2557

<xarray.Dataset>
Dimensions:          (RF: 136, YG: 400, XG: 1420, YC: 400, XC: 1420, time: 92, 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] 2005-10-01 2005-10-02 ... 2005-12-31
  * RC               (RC) float64 -1.25 -3.75 -6.25 ... -824.4 -881.7 -944.4
Data variables: (12/13)
    theta            (time, RC, YC, XC) float32 dask.array<chunksize=(1, 136, 400, 500), meta=np.ndarray>
    u                (time, RC, YC, XG) float32 dask.array<chunksize=(1, 136, 400, 500), meta=np.ndarray>
    v                (time, RC, YG, XC) float32 dask.array<chunksize=(1, 136, 400, 500), meta=np.ndarray>
    w                (time, RF, YC, XC) float32 dask.array<chunksize=(1, 136, 400, 500), meta=np.ndarray>
    salt             (time, RC, YC, XC) float32 dask.array<chunksize=(1, 136, 400, 500), meta=np.ndarray>
    KPP_diffusivity  (time, RF, YC, XC) float32 dask.array<chunksize=(1, 136, 400, 500), meta=np.ndarray>
    ...               ...
    N2               (time, RF, YC, XC) float32 dask.array<chunksize=(1, 1, 400, 500), meta=np.ndarray>
    S2               (time, RF, YC, XC) float32 dask.array<chunksize=(1, 1, 399, 499), meta=np.ndarray>
    Ri               (time, RF, YC, XC) float32 dask.array<chunksize=(1, 1, 399, 499), meta=np.ndarray>
    mld              (time, YC, XC) float64 dask.array<chunksize=(1, 400, 500), meta=np.ndarray>
    dcl_base         (time, YC, XC) float64 dask.array<chunksize=(1, 399, 499), meta=np.ndarray>
    dcl              (time, YC, XC) float64 dask.array<chunksize=(1, 399, 499), meta=np.ndarray>
Attributes:
    title:              daily snapshot from 1/20 degree Equatorial Pacific MI...
    easting:            longitude
    northing:           latitude
    field_julian_date:  488688
    julian_day_unit:    days since 1950-01-01 00:00:00
Hide code cell content 
ond08, metrics, grid = pump.model.read_mitgcm_20_year(
    state=True,
    start="2008-10-01",
    stop="2008-12-31",  # chunks={"RC": 10, "RF": 10},
)
ond08

3562 3653

<xarray.Dataset>
Dimensions:          (RF: 136, YG: 400, XG: 1420, YC: 400, XC: 1420, time: 92, 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-10-01 2008-10-02 ... 2008-12-31
  * RC               (RC) float64 -1.25 -3.75 -6.25 ... -824.4 -881.7 -944.4
Data variables: (12/13)
    theta            (time, RC, YC, XC) float32 dask.array<chunksize=(1, 136, 400, 500), meta=np.ndarray>
    u                (time, RC, YC, XG) float32 dask.array<chunksize=(1, 136, 400, 500), meta=np.ndarray>
    v                (time, RC, YG, XC) float32 dask.array<chunksize=(1, 136, 400, 500), meta=np.ndarray>
    w                (time, RF, YC, XC) float32 dask.array<chunksize=(1, 136, 400, 500), meta=np.ndarray>
    salt             (time, RC, YC, XC) float32 dask.array<chunksize=(1, 136, 400, 500), meta=np.ndarray>
    KPP_diffusivity  (time, RF, YC, XC) float32 dask.array<chunksize=(1, 136, 400, 500), meta=np.ndarray>
    ...               ...
    N2               (time, RF, YC, XC) float32 dask.array<chunksize=(1, 1, 400, 500), meta=np.ndarray>
    S2               (time, RF, YC, XC) float32 dask.array<chunksize=(1, 1, 399, 499), meta=np.ndarray>
    Ri               (time, RF, YC, XC) float32 dask.array<chunksize=(1, 1, 399, 499), meta=np.ndarray>
    mld              (time, YC, XC) float64 dask.array<chunksize=(1, 400, 500), meta=np.ndarray>
    dcl_base         (time, YC, XC) float64 dask.array<chunksize=(1, 399, 499), meta=np.ndarray>
    dcl              (time, YC, XC) float64 dask.array<chunksize=(1, 399, 499), meta=np.ndarray>
Attributes:
    title:              daily snapshot from 1/20 degree Equatorial Pacific MI...
    easting:            longitude
    northing:           latitude
    field_julian_date:  514992
    julian_day_unit:    days since 1950-01-01 00:00:00
Hide code cell content 
ond0508 = xr.concat([ond05.drop("time"), ond08.drop("time")], dim="year")
ond0508["year"] = [2005, 2008]

u, v, T#

I am looking for what’s different in 2005 vs 2008.

  1. EUC is in basically the same place. similar speed. similar latitude for EUC maximum.

  2. SEC is narrower in 2005 and faster.

  3. Cold tongue using top 100m mean temp is in basically the same place

  4. In 2005, \(v^2\) has a larger northern maximum. This is the largest difference I’ve seen so far

  5. MLD looks similar; deepest is 40m at eq

I think I want a metric for “TIW centroid” about the equator to judge if TIWs are further north during one season vs the other.

ds140 = dcpy.dask.map_copy(
    ond0508[["theta", "u", "v", "w"]].cf.sel(X=-140, method="nearest")
).load()

fg = (
    ds140.u.cf.sel(Z=slice(-250), YC=slice(-6, 6))
    .mean("time")
    .plot.contourf(col="year", robust=True, levels=np.arange(-1, 1.01, 0.05))
)

fg.map(lambda: plt.axvline(0, color="w"))

<xarray.plot.facetgrid.FacetGrid at 0x2af893d51640>
_images/9acb553a7d8cbc36036d28f2bc9cabbeff4cf7b18f98773154ee855680eab1c6.png 
argmax = ds140.u.sel(YC=slice(-2.5, 2.5)).argmax(["YC", "RC"])
ds140.u.sel(YC=slice(-2.5, 2.5)).YC.isel(YC=argmax["YC"]).plot(
    hue="year", size=5, aspect=2.5
)
plt.ylabel("Lat of EUC max")

Text(0, 0.5, 'Lat of EUC max')
_images/32729fd7f75a66216cb13c4c47590ed995d8fcfe281f80c6271f0b53bc2b04f7.png 
f, ax = plt.subplots(1, 2, sharey=True)

ds140.theta.cf.sel(Z=slice(-50)).cf.mean("Z").mean("time").plot(
    hue="year", y="YC", ax=ax[0]
)
(ds140.v**2 / 2).cf.sel(Z=slice(-50)).cf.mean(["Z", "time"]).plot(hue="year", y="YG")
ax[1].set_xlabel("$v^2/2$")

Text(0.5, 0, '$v^2/2$')
_images/d7402cf11e901a2f415efa14064ce005633201db6f4fa06074b3415c420fa9c5.png 
ds140.theta.cf.sel(Z=slice(-50)).cf.mean("Z").plot(
    robust=True, row="year", aspect=2, x="time"
)

<xarray.plot.facetgrid.FacetGrid at 0x2af82cf3a580>
_images/ac7241eadf981d3950c89fb7d67ccfba00bd7d9e106055c51051846d5f13811a.png

Cold tongue is in basically the same place

MLD#

meanmld = (
    ond0508.mld.sel(XC=-140, method="nearest")
    .sel(YC=slice(-0.5, 0.5))
    .mean("YC")
    .load()
)

meanmld.plot(hue="year", size=5, aspect=2.5)
plt.title("140°W, 0.5°S-0.5°N mean")

Text(0.5, 1.0, '140°W, 0.5°S-0.5°N mean')
_images/e674039877c3478b154c50151e6e106fd293a4146b9151170ba75a84905f566e.png

w#

Note near-equatorial downwelling in 2008

ds140.w.sel(RF=slice(-250), YC=slice(-6, 6)).mean("time").plot(col="year", robust=True)

<xarray.plot.facetgrid.FacetGrid at 0x2afdda761df0>
_images/56e3a04796baf8326f4e9b05ea735d53db8b1bf8a43394433ddfa4c9d51a5c9f.png 
def plot_w(w):
    slider = pnw.DiscreteSlider  # pnw.IntSlider(start=0, end=92)

    return (
        w.sel(RF=-50, method="nearest")
        .interactive.sel(time=slider)
        .hvplot(clim=(-3e-4, 3e-4), cmap="RdBu_r")
    )


wplot08 = plot_w(ond08.w)
wplot05 = plot_w(ond05.w)

wplot05 * (hv.HLine(0) * hv.VLine(-140))

wplot08 * (hv.HLine(0) * hv.VLine(-140))

w05 = ond05.w.sel(XC=[-155, -140, -125, -110], method="nearest").load()
w08 = ond08.w.sel(XC=[-155, -140, -125, -110], method="nearest").load()

w histograms by latitude#

  1. It looked like 2005 had more positive w; BUT once I change to log scale on the y-axis; it looks like 2008 has a lot more intense -ve w.

  2. The more intense -ve w is in 0.5S - 0.5N

  3. Are the TIWs in different places?

    1. At first I thought regridding to a y_EUC centered coordinate would help but it looks like the EUC max is in roughly the same place.

    2. I think I need a measure of TIW centroids or something?

Hide code cell source 
def plot_w_hist(*w, lon):
    kwargs = dict(bins=np.linspace(-4e-4, 4e-4, 100), histtype="step", density=True)
    breaks = [-3, -2, -1, -0.5, -0.25, 0.25, 0.5, 1, 2, 3]

    f, axes = plt.subplots(3, 3, sharex=True, sharey=True, constrained_layout=True)

    for ax, (left, right) in zip(axes.flat, zip(breaks[:-1], breaks[1:])):
        sel = lambda x: x.cf.sel(X=lon, method="nearest").sel(
            YC=slice(left, right), RF=slice(-100)
        )

        for w_ in w:
            w_.pipe(sel).plot.hist(
                **kwargs,
                label=str(w_.time.dt.year.data[0]),
                ax=ax,
                yscale="log",
                ylim=(1e1, 1e4),
            )
        ax.set_title(f"[{left!r}, {right!r})")

    f.legend(bbox_to_anchor=(1.15, 0.8))
    f.suptitle(f"lon={lon}")
    dcpy.plots.clean_axes(ax)


plot_w_hist(w05, w08, lon=-140)
_images/d0a3005f1072b5a23d2965adfb0cd73e39144dc36eadccc35160c83ef77935cf.png

TIWs, divergence terms, w#

  1. The w field is messy

  2. Smoothed v and ∂v/∂y line up nicely

w05.cf.sel(Z=slice(-50, -100)).sel(XC=-140, method="nearest").cf.mean("Z").plot(
    x="time", robust=True, size=4, aspect=3
)
w08.cf.sel(Z=slice(-50, -100)).sel(XC=-140, method="nearest").cf.mean("Z").plot(
    x="time", robust=True, size=4, aspect=3
)

<matplotlib.collections.QuadMesh at 0x2afd9360b5e0>
_images/b2cd50604c7d97a58f4ff11225a5ea53facb48ba6825e9442d4947992450c6d8.png _images/b800f90e82f155aece89235787cc68b59012c7781e8780459cd419579197f5d5.png

v#

shows nice phase signature and TIWs “leaning” into the SEC-EUC zonal flow

fg = (
    ds140.v.cf.sel(Z=slice(-30, -100))
    .cf.mean("Z")
    .rolling(time=5, center=True)
    .mean()
    .plot(robust=True, row="year", aspect=2, x="time", cmap=mpl.cm.coolwarm)
)
_images/f3ce2815bd566dc1aace50563ab796c6a33330937459dad4556aef231b494c1f.png

∂v/∂y#

fg = (
    ds140.v.cf.differentiate("Y")
    .cf.sel(Z=slice(-30, -100))
    .cf.mean("Z")
    .rolling(time=5, center=True)
    .mean()
    .plot(robust=True, row="year", aspect=2, x="time", cmap=mpl.cm.coolwarm)
)
_images/8f591c1a06cdfcf5fd30af9892a1d6aab3b72fafb9fe6dcf52e80851880a14c2.png

w#

smoothing the w field helps a bit; see downwelling followed by upwelling TIW pattern

fg = (
    ds140.w.cf.sel(Z=slice(-40, -120))
    .cf.mean("Z")
    .rolling(time=3, YC=3, center=True)
    .mean()
    .plot(robust=True, row="year", aspect=2, x="time", cmap=mpl.cm.coolwarm)
)
fg.data = (
    ds140.v.cf.differentiate("Y")
    .cf.sel(Z=slice(-40, -100))
    .cf.mean("Z")
    .rolling(time=5, YG=5, center=True)
    .mean()
)
# fg.map_dataarray(
#    xr.plot.contour,
#    x="time",
#    y="YG",
#    levels=5,
#    cmap=mpl.cm.PuOr_r,
#    linewidths=1,
#    robust=True,
#    add_colorbar=False,
# )
_images/b642dbe5a967438d6a6e4f70f709216bdbfde2c8a9c2848fed786c9068ac65df.png

I started looking to see if there was a dominant contributor to \(u_x + v_y\) but not clear

ux = grid.derivative(ond05.isel(time=50).u, "X").load()
vy = grid.derivative(ond05.isel(time=50).v, "Y").load()

/glade/u/home/dcherian/python/xgcm/xgcm/grid.py:1384: UserWarning: Metric at ('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 ('RC', 'YC', 'XC') being interpolated from metrics at dimensions ('YG', 'XC'). Boundary value set to 'extend'.
  warnings.warn(

div = xr.concat([ux, vy], dim="div")
div["div"] = ["ux", "vy"]

(ux + vy).sel(RC=-150, method="nearest").plot(
    robust=True,
    add_colorbar=True,
    # vmin=-1,
    # vmax=1,
    cmap=mpl.cm.RdBu_r,
    aspect=2,
    size=4,
    ylim=(-2.5, 6),
)

<matplotlib.collections.QuadMesh at 0x2afd91decd90>
_images/14f5032875f18d131be18ef8790798623f5555cf0dc3831b9eb2c35566583837.png

Single TIW period: 2021 paper#

period4 = xr.open_zarr(
    "/glade/work/dcherian/pump/zarrs/110w-period-4-3.zarr", consolidated=True
)
del period4["time"].attrs["long_name"]
central = period4.isel(longitude=1)

metrics = (
    pump.model.read_metrics(
        f"/glade/campaign/cgd/oce/people/dcherian/TPOS_MITgcm_1_hb/"
    )
    .reindex_like(central.expand_dims("longitude"), method="nearest")
    .load()
).squeeze()

/glade/u/home/dcherian/miniconda3/envs/dcpy/lib/python3.8/site-packages/xarray/core/indexing.py:1375: 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,)]

central["Jq"] = (
    1035 * 3994 * (central.DFrI_TH + central.KPPg_TH.fillna(0)) / metrics.RAC
)
central["time"] = central["time"] - pd.Timedelta("7h")
central["dens"] = pump.mdjwf.dens(central.salt, central.theta, np.array([0.0]))
central = pump.calc.calc_reduced_shear(central)
central["mld"] = pump.calc.get_mld(central.dens)
central["dcl_base"] = pump.calc.get_dcl_base_Ri(central)
central["dcl"] = central["mld"] - central["dcl_base"]
central["sst"] = central.theta.isel(depth=0)
central["eucmax"] = pump.calc.get_euc_max(
    central.u.sel(latitude=0, method="nearest", drop=True)
)
central.depth.attrs["units"] = "m"
central.Jq.attrs["long_name"] = "$J^t_q$"
central.Jq.attrs["units"] = "W/m²"
central.S2.attrs["long_name"] = "$S²$"
central.S2.attrs["units"] = "$s^{-2}$"
central.Ri.attrs["long_name"] = "$Ri$"
central.N2.attrs["long_name"] = "$N²$"
central.N2.attrs["units"] = "$s^{-2}$"
central

calc uz
calc vz
calc S2
calc N2
calc shred2
Calc Ri

<xarray.Dataset>
Dimensions:    (depth: 200, latitude: 400, time: 779)
Coordinates:
  * depth      (depth) float32 -0.5 -1.5 -2.5 -3.5 ... -197.5 -198.5 -199.5
  * latitude   (latitude) float32 -10.0 -9.95 -9.9 -9.85 ... 9.85 9.9 9.95 10.0
    longitude  float32 -110.0
  * time       (time) datetime64[ns] 1995-11-14T17:00:00 ... 1995-12-17T03:00:00
Data variables:
    DFrI_TH    (time, depth, latitude) float32 dask.array<chunksize=(1, 200, 400), meta=np.ndarray>
    ETAN       (time, latitude) float32 dask.array<chunksize=(1, 400), meta=np.ndarray>
    KPPRi      (time, depth, latitude) float32 dask.array<chunksize=(1, 200, 400), meta=np.ndarray>
    KPPbo      (time, latitude) float32 dask.array<chunksize=(1, 400), meta=np.ndarray>
    KPPdiffT   (time, depth, latitude) float32 dask.array<chunksize=(1, 200, 400), meta=np.ndarray>
    KPPg_TH    (time, depth, latitude) float32 dask.array<chunksize=(1, 200, 400), meta=np.ndarray>
    KPPhbl     (time, latitude) float32 dask.array<chunksize=(1, 400), meta=np.ndarray>
    KPPviscA   (time, depth, latitude) float32 dask.array<chunksize=(1, 200, 400), meta=np.ndarray>
    Um_Diss    (time, depth, latitude) float32 dask.array<chunksize=(1, 200, 400), meta=np.ndarray>
    VISrI_Um   (time, depth, latitude) float32 dask.array<chunksize=(1, 200, 400), meta=np.ndarray>
    VISrI_Vm   (time, depth, latitude) float32 dask.array<chunksize=(1, 200, 400), meta=np.ndarray>
    Vm_Diss    (time, depth, latitude) float32 dask.array<chunksize=(1, 200, 400), meta=np.ndarray>
    salt       (time, depth, latitude) float32 dask.array<chunksize=(1, 200, 400), meta=np.ndarray>
    theta      (time, depth, latitude) float32 dask.array<chunksize=(1, 200, 400), meta=np.ndarray>
    u          (time, depth, latitude) float32 dask.array<chunksize=(1, 200, 400), meta=np.ndarray>
    v          (time, depth, latitude) float32 dask.array<chunksize=(1, 200, 400), meta=np.ndarray>
    w          (time, depth, latitude) float32 dask.array<chunksize=(1, 200, 400), meta=np.ndarray>
    Jq         (time, depth, latitude) float64 dask.array<chunksize=(1, 200, 400), meta=np.ndarray>
    dens       (time, depth, latitude) float64 dask.array<chunksize=(1, 200, 400), meta=np.ndarray>
    uz         (time, depth, latitude) float32 dask.array<chunksize=(1, 200, 400), meta=np.ndarray>
    vz         (time, depth, latitude) float32 dask.array<chunksize=(1, 200, 400), meta=np.ndarray>
    S2         (time, depth, latitude) float32 dask.array<chunksize=(1, 200, 400), meta=np.ndarray>
    shear      (time, depth, latitude) float32 dask.array<chunksize=(1, 200, 400), meta=np.ndarray>
    N2         (time, depth, latitude) float64 dask.array<chunksize=(1, 200, 400), meta=np.ndarray>
    shred2     (time, depth, latitude) float64 dask.array<chunksize=(1, 200, 400), meta=np.ndarray>
    Ri         (time, depth, latitude) float64 dask.array<chunksize=(1, 200, 400), meta=np.ndarray>
    mld        (time, latitude) float32 dask.array<chunksize=(1, 400), meta=np.ndarray>
    dcl_base   (time, latitude) float32 dask.array<chunksize=(1, 400), meta=np.ndarray>
    dcl        (time, latitude) float32 dask.array<chunksize=(1, 400), meta=np.ndarray>
    sst        (time, latitude) float32 dask.array<chunksize=(1, 400), meta=np.ndarray>
    eucmax     (time) float32 -97.5 -97.5 -97.5 -97.5 ... -100.5 -99.5 -99.5
Attributes:
    easting:            longitude
    field_julian_date:  400296
    julian_day_unit:    days since 1950-01-01 00:00:00
    northing:           latitude
    title:              daily snapshot from 1/20 degree Equatorial Pacific MI...
central.w.load()

<xarray.DataArray 'w' (time: 779, depth: 200, latitude: 400, longitude: 3)>
array([[[[-8.32392686e-07, -8.61631122e-07, -9.07944013e-07],
         [-8.17283251e-07, -8.39938480e-07, -8.78521575e-07],
         [-8.14519808e-07, -8.35834385e-07, -8.55945700e-07],
         ...,
         [-2.07602127e-07, -2.00864179e-07, -1.98754080e-07],
         [-2.06921740e-07, -2.26620955e-07, -2.43814696e-07],
         [-1.80513695e-07, -2.17974403e-07, -2.52572022e-07]],

        [[-2.42294459e-06, -3.05959338e-06, -3.24233270e-06],
         [-2.51864435e-06, -3.08625977e-06, -3.58562716e-06],
         [-2.49380400e-06, -3.41150280e-06, -3.97996473e-06],
         ...,
         [ 2.10747203e-06,  1.65012580e-06,  1.15887929e-06],
         [ 1.27279213e-06,  7.15527619e-07,  2.72692631e-07],
         [-5.73148554e-07, -1.04561616e-06, -1.42848603e-06]],

        [[-4.02487240e-06, -5.26639315e-06, -5.58513466e-06],
         [-4.23122447e-06, -5.34199489e-06, -6.30122440e-06],
         [-4.18542413e-06, -5.99460918e-06, -7.11171651e-06],
         ...,
...
         [ 4.57714168e-06,  3.17592526e-06, -1.33422600e-05],
         [-1.43890775e-05, -3.12485827e-05, -5.46250267e-05],
         [-5.27998491e-05, -7.76431771e-05, -1.01270351e-04]],

        [[ 3.85325366e-05,  2.63559486e-05,  1.21922913e-05],
         [ 2.41746147e-05,  1.80613529e-06, -2.40914906e-05],
         [ 5.65520440e-06, -2.40457211e-05, -4.94654341e-05],
         ...,
         [ 4.26078304e-06,  2.97869406e-06, -1.33218573e-05],
         [-1.47053061e-05, -3.17766717e-05, -5.48467360e-05],
         [-5.34974752e-05, -7.84654549e-05, -1.01884856e-04]],

        [[ 4.04056191e-05,  2.77295694e-05,  1.29975933e-05],
         [ 2.59654771e-05,  3.19075821e-06, -2.32989441e-05],
         [ 6.97214091e-06, -2.29840607e-05, -4.85406090e-05],
         ...,
         [ 3.94934204e-06,  2.77228696e-06, -1.33081094e-05],
         [-1.50320866e-05, -3.23151435e-05, -5.50692130e-05],
         [-5.41884583e-05, -7.92831343e-05, -1.02494916e-04]]]],
      dtype=float32)
Coordinates:
  * depth      (depth) float32 -0.5 -1.5 -2.5 -3.5 ... -197.5 -198.5 -199.5
  * latitude   (latitude) float32 -10.0 -9.95 -9.9 -9.85 ... 9.85 9.9 9.95 10.0
  * longitude  (longitude) float32 -110.1 -110.0 -110.0
  * time       (time) datetime64[ns] 1995-11-15 ... 1995-12-17T10:00:00

Maps#

lat-lon at depth#

(
    central.w.sel(depth=np.arange(-5, -100, -20), method="nearest").plot(
        robust=True, x="time", col="depth", col_wrap=3
    )
)

<xarray.plot.facetgrid.FacetGrid at 0x2b12cc356310>
_images/f976bcac0c2e6005b3c1407e11c94628b7758cf34df6a246e2fb1962bb218b71.png

lat-lon top 85m mean#

(central.w.sel(depth=slice(-84.5)).mean("depth").plot(robust=True, x="time"))

<matplotlib.collections.QuadMesh at 0x2b12cc585730>
_images/54872aa39e5925b9b836a978928c3d4317f27b2e5d8123a5fc5ccb42f38d9dff.png

depth-time at lats#

lats = [4, 3, 2, 1, 0, -1]

central.theta.isel(depth=0).plot(robust=True, x="time", cmap=mpl.cm.RdYlBu_r)
[plt.axhline(lat) for lat in lats]

(
    central.w.sel(latitude=lats, method="nearest").plot(
        col_wrap=2,
        col="latitude",
        x="time",
        robust=True,
        cbar_kwargs={"orientation": "horizontal"},
    )
)

<xarray.plot.facetgrid.FacetGrid at 0x2b12cea61610>
_images/8f096f9b2fa7372f39d07338af41070987403d010f2371e57a9840ab53da487e.png _images/47cfca5e7b838a4f8ec66451437d4ec68eba1d2a1346d9563f01f1a29c73b482.png

depth-lat at times#

times = [
    "1995-11-17",
    "1995-11-25",
    "1995-11-26",
    "1995-11-27",
    "1995-12-01",
    "1995-12-10",
]

central.theta.isel(depth=0).plot(robust=True, x="time", cmap=mpl.cm.RdYlBu_r)
[plt.axvline(time) for time in times]

fg = central.w.sel(time=[pd.Timestamp(tt) for tt in times]).plot(
    y="depth",
    row="time",
    robust=True,
    cbar_kwargs={"orientation": "horizontal"},
    xlim=(-5, 5),
)
fg.fig.set_size_inches((8, 6))
_images/b0a038c49983b5a1fe3b3a2948f7237632c18bdeed3d0507239f6c045fff7e9c.png _images/4b801a4e09d56778a836e5ccf6569e5609d8a5a833b400efb83800a024f3e972.png

Spectra#

import dcpy

Why is the spectrum white for freq < 2 days?

z = -50
lats = [0, 3, 3.5, 4]

subset = central.w.reindex(latitude=lats, method="nearest").sel(
    depth=z, method="nearest"
)

for lat in lats:
    dcpy.ts.PlotSpectrum(subset.sel(latitude=lat), label=str(lat), ax=plt.gca())

plt.legend(title="Lat (°N)")
plt.title(f"depth = {z}m")

Text(0.5, 1.0, 'depth = -50m')
_images/fbce46f5c0de60f1bf1a5dc7049bf52c01b2adeb80a69322cc130e728def4df0.png 
z = -150

subset = central.w.reindex(latitude=lats, method="nearest").sel(
    depth=z, method="nearest"
)

for lat in lats:
    dcpy.ts.PlotSpectrum(subset.sel(latitude=lat), label=str(lat), ax=plt.gca())

plt.legend(title="Lat (°N)")
plt.title(f"depth = {z}m")

Text(0.5, 1.0, 'depth = -150m')
_images/e34a1737cc6ee9934089a97b2c09c049c521f90b7036e49d49ae4c89e85b7154.png

Mixing in T space#

central = central.persist()

plt.rcParams["figure.figsize"] = (8.0, 5.0)
plt.rcParams["figure.dpi"] = 140

def plot_Jq_T(subset):
    (-1 * subset.Jq.rolling(depth=2).mean()).compute().differentiate("depth").plot(
        x="time", vmin=-100, vmax=100, cmap=mpl.cm.RdBu_r
    )
    subset.theta.plot.contour(x="time", levels=[20, 21, 22, 23, 23.5, 24])
    subset.eucmax.plot.line(x="time", color="k", ls="--")
    plt.gca().set_ylim((-125, 0))
    plt.gcf().set_size_inches((8, 4))


plot_Jq_T(central.sel(latitude=3.5, method="nearest"))
plt.figure()
plot_Jq_T(central.sel(latitude=0, method="nearest"))
_images/4fdfc5092be7bab006c68ced116233c4033c4adbbf84d5a1895503d832d7b527.png _images/d6886582afc402b7b2f7302839fd8e76649ca8c504399766b4bac291d2722f78.png 
client.cancel(central)