Reading ECCO \(χ_{redi}\)#

%load_ext watermark

from pathlib import Path

import cf_xarray as cfxr
import dask
import ecco_v4_py as ecco
import matplotlib as mpl
import ncar_jobqueue
import numpy as np
import xgcm
from distributed import Client

import xarray as xr

%watermark -iv

xr.set_options(display_expand_data=False)

dirname = "/glade/work/dcherian/mitgcm/ECCOV4/release4/run/diags/"


def open_directory(dirname):
    import xmitgcm

    ds = xr.merge(
        [
            dask.optimize(
                xmitgcm.open_mdsdataset(
                    str(path),
                    dirname + "../",
                    geometry="llc",
                    chunks={"face": 1, "k_l": -1, "k": -1},
                    llc_method="smallchunks",
                    read_grid=idx == 0,
                    default_dtype=np.float32,
                )
                .isel(time=slice(-12))
                .chunk({"time": 120})
            )[0]
            for idx, path in enumerate(Path(dirname).glob("*"))
        ],
        compat="override",
    )  # .sel(face=2)
    ds.coords["drC_kl"] = ds.drC.isel(k_p1=slice(-1)).rename({"k_p1": "k_l"}).variable

    grid = xgcm.Grid(
        ds,
        periodic=False,
        boundary="fill",
        fill_value=np.nan,
        metrics={
            ("X",): ["dxC", "dxG"],  # X distances
            ("Y",): ["dyC", "dyG"],  # Y distances
            ("Z",): ["drF", "drC", "drC_kl"],  # Z distances
            ("X", "Y"): ["rA", "rAz", "rAs", "rAw"],  # Areas
        },
    )
    ds["Tx"] = grid.derivative(ds.THETA, "X")
    ds["Ty"] = grid.derivative(ds.THETA, "Y")
    ds["Tz"] = -1 * grid.derivative(ds.THETA, "Z")

    ds["GM_CHI"] = (
        grid.interp(ds.GM_CHIX + ds.GM_CHIXO, "X", to="center")
        + grid.interp(ds.GM_CHIY + ds.GM_CHIYO, "Y", to="center")
        + grid.interp(ds.GM_CHIZ, "Z", to="center")
    )
    return ds, grid

The watermark extension is already loaded. To reload it, use:
  %reload_ext watermark
pandas       : 1.5.3
cf_xarray    : 0.8.0
dask         : 2023.3.2
numpy        : 1.23.5
ncar_jobqueue: 2021.4.14
xgcm         : 0.6.1
json         : 2.0.9
ecco_v4_py   : 1.5.5
xmitgcm      : 0.5.2
sys          : 3.10.10 | packaged by conda-forge | (main, Mar 24 2023, 20:08:06) [GCC 11.3.0]
matplotlib   : 3.7.1
xarray       : 2023.3.0

cluster = ncar_jobqueue.NCARCluster(
    local_directory="/local_scratch/pbs.$PBS_JOBID/dask/spill"
)
cluster
cluster.adapt(minimum_jobs=1, maximum_jobs=32)

/glade/u/home/dcherian/miniconda3/envs/eddydiff/lib/python3.10/site-packages/dask_jobqueue/core.py:255: FutureWarning: job_extra has been renamed to job_extra_directives. You are still using it (even if only set to []; please also check config files). If you did not set job_extra_directives yet, job_extra will be respected for now, but it will be removed in a future release. If you already set job_extra_directives, job_extra is ignored and you can remove it.
  warnings.warn(warn, FutureWarning)
/glade/u/home/dcherian/miniconda3/envs/eddydiff/lib/python3.10/site-packages/dask_jobqueue/core.py:274: FutureWarning: env_extra has been renamed to job_script_prologue. You are still using it (even if only set to []; please also check config files). If you did not set job_script_prologue yet, env_extra will be respected for now, but it will be removed in a future release. If you already set job_script_prologue, env_extra is ignored and you can remove it.
  warnings.warn(warn, FutureWarning)
/glade/u/home/dcherian/miniconda3/envs/eddydiff/lib/python3.10/site-packages/dask_jobqueue/core.py:255: FutureWarning: job_extra has been renamed to job_extra_directives. You are still using it (even if only set to []; please also check config files). If you did not set job_extra_directives yet, job_extra will be respected for now, but it will be removed in a future release. If you already set job_extra_directives, job_extra is ignored and you can remove it.
  warnings.warn(warn, FutureWarning)
/glade/u/home/dcherian/miniconda3/envs/eddydiff/lib/python3.10/site-packages/dask_jobqueue/core.py:274: FutureWarning: env_extra has been renamed to job_script_prologue. You are still using it (even if only set to []; please also check config files). If you did not set job_script_prologue yet, env_extra will be respected for now, but it will be removed in a future release. If you already set job_script_prologue, env_extra is ignored and you can remove it.
  warnings.warn(warn, FutureWarning)

<distributed.deploy.adaptive.Adaptive at 0x2abdfae84310>

client = Client(cluster)
client

Client

Client-56c116a5-e061-11ed-bfe1-3cecef1acc54

Connection method: Cluster object Cluster type: dask_jobqueue.PBSCluster
Dashboard: https://jupyterhub.hpc.ucar.edu/stable/user/dcherian/proxy/8787/status

Cluster Info

 
%ls /glade/work/dcherian/mitgcm/ECCOV4/release4/run/diags/

GM_CHIX/   GM_CHIY/   GM_CHIZ/  GM_Kuz/  GM_Kvz/  GM_Kwy/  GM_KwzTz/  THETA/
GM_CHIXO/  GM_CHIYO/  GM_Kux/   GM_Kvy/  GM_Kwx/  GM_Kwz/  GM_ubT/

Compare#

adj, grid = open_directory("/glade/work/dcherian/mitgcm/ECCOV4/release4/run/diags/")
const, _ = open_directory("/glade/work/dcherian/mitgcm/ECCOV4/release4/run2/diags/")
display(adj)
display(const)

<xarray.Dataset>
Dimensions:    (i: 90, i_g: 90, j: 90, j_g: 90, k: 50, k_u: 50, k_l: 50,
                k_p1: 51, face: 13, time: 176)
Coordinates: (12/45)
  * i          (i) int64 0 1 2 3 4 5 6 7 8 9 ... 80 81 82 83 84 85 86 87 88 89
  * i_g        (i_g) int64 0 1 2 3 4 5 6 7 8 9 ... 80 81 82 83 84 85 86 87 88 89
  * j          (j) int64 0 1 2 3 4 5 6 7 8 9 ... 80 81 82 83 84 85 86 87 88 89
  * j_g        (j_g) int64 0 1 2 3 4 5 6 7 8 9 ... 80 81 82 83 84 85 86 87 88 89
  * k          (k) int64 0 1 2 3 4 5 6 7 8 9 ... 40 41 42 43 44 45 46 47 48 49
  * k_u        (k_u) int64 0 1 2 3 4 5 6 7 8 9 ... 40 41 42 43 44 45 46 47 48 49
    ...         ...
    maskCtrlC  (k, face, j, i) bool dask.array<chunksize=(50, 1, 90, 90), meta=np.ndarray>
    maskCtrlS  (k, face, j_g, i) bool dask.array<chunksize=(50, 1, 90, 90), meta=np.ndarray>
    rhoRef     (k) >f4 dask.array<chunksize=(50,), meta=np.ndarray>
    maskCtrlW  (k, face, j, i_g) bool dask.array<chunksize=(50, 1, 90, 90), meta=np.ndarray>
    iter       (time) float64 dask.array<chunksize=(120,), meta=np.ndarray>
    drC_kl     (k_l) >f4 dask.array<chunksize=(50,), meta=np.ndarray>
Data variables: (12/19)
    GM_CHIXO   (time, k, face, j, i_g) float32 dask.array<chunksize=(120, 50, 1, 90, 90), meta=np.ndarray>
    GM_CHIZ    (time, k_l, face, j, i) float32 dask.array<chunksize=(120, 50, 1, 90, 90), meta=np.ndarray>
    GM_ubT     (time, k, face, j, i_g) float32 dask.array<chunksize=(120, 50, 1, 90, 90), meta=np.ndarray>
    GM_CHIX    (time, k, face, j, i_g) float32 dask.array<chunksize=(120, 50, 1, 90, 90), meta=np.ndarray>
    GM_CHIYO   (time, k, face, j_g, i) float32 dask.array<chunksize=(120, 50, 1, 90, 90), meta=np.ndarray>
    GM_Kvz     (time, k, face, j_g, i) float32 dask.array<chunksize=(120, 50, 1, 90, 90), meta=np.ndarray>
    ...         ...
    GM_CHIY    (time, k, face, j_g, i) float32 dask.array<chunksize=(120, 50, 1, 90, 90), meta=np.ndarray>
    GM_Kuz     (time, k, face, j, i_g) float32 dask.array<chunksize=(120, 50, 1, 90, 90), meta=np.ndarray>
    Tx         (time, k, face, j, i_g) float32 dask.array<chunksize=(120, 50, 1, 90, 1), meta=np.ndarray>
    Ty         (time, k, face, j_g, i) float32 dask.array<chunksize=(120, 50, 1, 1, 90), meta=np.ndarray>
    Tz         (time, k_l, face, j, i) float32 dask.array<chunksize=(120, 1, 1, 90, 90), meta=np.ndarray>
    GM_CHI     (time, k, face, j, i) float32 dask.array<chunksize=(120, 49, 1, 89, 89), meta=np.ndarray>
Attributes:
    Conventions:  CF-1.6
    title:        netCDF wrapper of MITgcm MDS binary data
    source:       MITgcm
    history:      Created by calling `open_mdsdataset(grid_dir='/glade/work/d...
<xarray.Dataset>
Dimensions:    (i: 90, i_g: 90, j: 90, j_g: 90, k: 50, k_u: 50, k_l: 50,
                k_p1: 51, face: 13, time: 176)
Coordinates: (12/45)
  * i          (i) int64 0 1 2 3 4 5 6 7 8 9 ... 80 81 82 83 84 85 86 87 88 89
  * i_g        (i_g) int64 0 1 2 3 4 5 6 7 8 9 ... 80 81 82 83 84 85 86 87 88 89
  * j          (j) int64 0 1 2 3 4 5 6 7 8 9 ... 80 81 82 83 84 85 86 87 88 89
  * j_g        (j_g) int64 0 1 2 3 4 5 6 7 8 9 ... 80 81 82 83 84 85 86 87 88 89
  * k          (k) int64 0 1 2 3 4 5 6 7 8 9 ... 40 41 42 43 44 45 46 47 48 49
  * k_u        (k_u) int64 0 1 2 3 4 5 6 7 8 9 ... 40 41 42 43 44 45 46 47 48 49
    ...         ...
    maskCtrlC  (k, face, j, i) bool dask.array<chunksize=(50, 1, 90, 90), meta=np.ndarray>
    maskCtrlS  (k, face, j_g, i) bool dask.array<chunksize=(50, 1, 90, 90), meta=np.ndarray>
    rhoRef     (k) >f4 dask.array<chunksize=(50,), meta=np.ndarray>
    maskCtrlW  (k, face, j, i_g) bool dask.array<chunksize=(50, 1, 90, 90), meta=np.ndarray>
    iter       (time) float64 dask.array<chunksize=(120,), meta=np.ndarray>
    drC_kl     (k_l) >f4 dask.array<chunksize=(50,), meta=np.ndarray>
Data variables: (12/19)
    GM_CHIXO   (time, k, face, j, i_g) float32 dask.array<chunksize=(120, 50, 1, 90, 90), meta=np.ndarray>
    GM_CHIZ    (time, k_l, face, j, i) float32 dask.array<chunksize=(120, 50, 1, 90, 90), meta=np.ndarray>
    GM_ubT     (time, k, face, j, i_g) float32 dask.array<chunksize=(120, 50, 1, 90, 90), meta=np.ndarray>
    GM_CHIX    (time, k, face, j, i_g) float32 dask.array<chunksize=(120, 50, 1, 90, 90), meta=np.ndarray>
    GM_CHIYO   (time, k, face, j_g, i) float32 dask.array<chunksize=(120, 50, 1, 90, 90), meta=np.ndarray>
    GM_Kvz     (time, k, face, j_g, i) float32 dask.array<chunksize=(120, 50, 1, 90, 90), meta=np.ndarray>
    ...         ...
    GM_CHIY    (time, k, face, j_g, i) float32 dask.array<chunksize=(120, 50, 1, 90, 90), meta=np.ndarray>
    GM_Kuz     (time, k, face, j, i_g) float32 dask.array<chunksize=(120, 50, 1, 90, 90), meta=np.ndarray>
    Tx         (time, k, face, j, i_g) float32 dask.array<chunksize=(120, 50, 1, 90, 1), meta=np.ndarray>
    Ty         (time, k, face, j_g, i) float32 dask.array<chunksize=(120, 50, 1, 1, 90), meta=np.ndarray>
    Tz         (time, k_l, face, j, i) float32 dask.array<chunksize=(120, 1, 1, 90, 90), meta=np.ndarray>
    GM_CHI     (time, k, face, j, i) float32 dask.array<chunksize=(120, 49, 1, 89, 89), meta=np.ndarray>
Attributes:
    Conventions:  CF-1.6
    title:        netCDF wrapper of MITgcm MDS binary data
    source:       MITgcm
    history:      Created by calling `open_mdsdataset(grid_dir='/glade/work/d...
concatted = xr.concat(
    [adj, const], dim="node", compat="override", coords="minimal"
).assign_coords(node=["adj", "const"])

<xarray.Dataset>
Dimensions:    (i: 90, i_g: 90, j: 90, j_g: 90, k: 50, k_u: 50, k_l: 50,
                k_p1: 51, face: 13, time: 176, node: 2)
Coordinates: (12/46)
  * i          (i) int64 0 1 2 3 4 5 6 7 8 9 ... 80 81 82 83 84 85 86 87 88 89
  * i_g        (i_g) int64 0 1 2 3 4 5 6 7 8 9 ... 80 81 82 83 84 85 86 87 88 89
  * j          (j) int64 0 1 2 3 4 5 6 7 8 9 ... 80 81 82 83 84 85 86 87 88 89
  * j_g        (j_g) int64 0 1 2 3 4 5 6 7 8 9 ... 80 81 82 83 84 85 86 87 88 89
  * k          (k) int64 0 1 2 3 4 5 6 7 8 9 ... 40 41 42 43 44 45 46 47 48 49
  * k_u        (k_u) int64 0 1 2 3 4 5 6 7 8 9 ... 40 41 42 43 44 45 46 47 48 49
    ...         ...
    maskCtrlS  (k, face, j_g, i) bool dask.array<chunksize=(50, 1, 90, 90), meta=np.ndarray>
    rhoRef     (k) >f4 dask.array<chunksize=(50,), meta=np.ndarray>
    maskCtrlW  (k, face, j, i_g) bool dask.array<chunksize=(50, 1, 90, 90), meta=np.ndarray>
    iter       (time) float64 dask.array<chunksize=(120,), meta=np.ndarray>
    drC_kl     (k_l) >f4 dask.array<chunksize=(50,), meta=np.ndarray>
  * node       (node) <U5 'adj' 'const'
Data variables: (12/19)
    GM_CHIXO   (node, time, k, face, j, i_g) float32 dask.array<chunksize=(1, 120, 50, 1, 90, 90), meta=np.ndarray>
    GM_CHIZ    (node, time, k_l, face, j, i) float32 dask.array<chunksize=(1, 120, 50, 1, 90, 90), meta=np.ndarray>
    GM_ubT     (node, time, k, face, j, i_g) float32 dask.array<chunksize=(1, 120, 50, 1, 90, 90), meta=np.ndarray>
    GM_CHIX    (node, time, k, face, j, i_g) float32 dask.array<chunksize=(1, 120, 50, 1, 90, 90), meta=np.ndarray>
    GM_CHIYO   (node, time, k, face, j_g, i) float32 dask.array<chunksize=(1, 120, 50, 1, 90, 90), meta=np.ndarray>
    GM_Kvz     (node, time, k, face, j_g, i) float32 dask.array<chunksize=(1, 120, 50, 1, 90, 90), meta=np.ndarray>
    ...         ...
    GM_CHIY    (node, time, k, face, j_g, i) float32 dask.array<chunksize=(1, 120, 50, 1, 90, 90), meta=np.ndarray>
    GM_Kuz     (node, time, k, face, j, i_g) float32 dask.array<chunksize=(1, 120, 50, 1, 90, 90), meta=np.ndarray>
    Tx         (node, time, k, face, j, i_g) float32 dask.array<chunksize=(1, 120, 50, 1, 90, 1), meta=np.ndarray>
    Ty         (node, time, k, face, j_g, i) float32 dask.array<chunksize=(1, 120, 50, 1, 1, 90), meta=np.ndarray>
    Tz         (node, time, k_l, face, j, i) float32 dask.array<chunksize=(1, 120, 1, 1, 90, 90), meta=np.ndarray>
    GM_CHI     (node, time, k, face, j, i) float32 dask.array<chunksize=(1, 120, 49, 1, 89, 89), meta=np.ndarray>
Attributes:
    Conventions:  CF-1.6
    title:        netCDF wrapper of MITgcm MDS binary data
    source:       MITgcm
    history:      Created by calling `open_mdsdataset(grid_dir='/glade/work/d...
natre = concatted.sel(face=2).isel(
    i=slice(10, 18), j=slice(10, 20), i_g=slice(10, 18), j_g=slice(10, 20)
)
natre

<xarray.Dataset>
Dimensions:    (i: 8, i_g: 8, j: 10, j_g: 10, k: 50, k_u: 50, k_l: 50,
                k_p1: 51, time: 176, node: 2)
Coordinates: (12/46)
  * i          (i) int64 10 11 12 13 14 15 16 17
  * i_g        (i_g) int64 10 11 12 13 14 15 16 17
  * j          (j) int64 10 11 12 13 14 15 16 17 18 19
  * j_g        (j_g) int64 10 11 12 13 14 15 16 17 18 19
  * k          (k) int64 0 1 2 3 4 5 6 7 8 9 ... 40 41 42 43 44 45 46 47 48 49
  * k_u        (k_u) int64 0 1 2 3 4 5 6 7 8 9 ... 40 41 42 43 44 45 46 47 48 49
    ...         ...
    maskCtrlS  (k, j_g, i) bool dask.array<chunksize=(50, 10, 8), meta=np.ndarray>
    rhoRef     (k) >f4 dask.array<chunksize=(50,), meta=np.ndarray>
    maskCtrlW  (k, j, i_g) bool dask.array<chunksize=(50, 10, 8), meta=np.ndarray>
    iter       (time) float64 dask.array<chunksize=(120,), meta=np.ndarray>
    drC_kl     (k_l) >f4 dask.array<chunksize=(50,), meta=np.ndarray>
  * node       (node) <U5 'adj' 'const'
Data variables: (12/19)
    GM_CHIXO   (node, time, k, j, i_g) float32 dask.array<chunksize=(1, 120, 50, 10, 8), meta=np.ndarray>
    GM_CHIZ    (node, time, k_l, j, i) float32 dask.array<chunksize=(1, 120, 50, 10, 8), meta=np.ndarray>
    GM_ubT     (node, time, k, j, i_g) float32 dask.array<chunksize=(1, 120, 50, 10, 8), meta=np.ndarray>
    GM_CHIX    (node, time, k, j, i_g) float32 dask.array<chunksize=(1, 120, 50, 10, 8), meta=np.ndarray>
    GM_CHIYO   (node, time, k, j_g, i) float32 dask.array<chunksize=(1, 120, 50, 10, 8), meta=np.ndarray>
    GM_Kvz     (node, time, k, j_g, i) float32 dask.array<chunksize=(1, 120, 50, 10, 8), meta=np.ndarray>
    ...         ...
    GM_CHIY    (node, time, k, j_g, i) float32 dask.array<chunksize=(1, 120, 50, 10, 8), meta=np.ndarray>
    GM_Kuz     (node, time, k, j, i_g) float32 dask.array<chunksize=(1, 120, 50, 10, 8), meta=np.ndarray>
    Tx         (node, time, k, j, i_g) float32 dask.array<chunksize=(1, 120, 50, 10, 8), meta=np.ndarray>
    Ty         (node, time, k, j_g, i) float32 dask.array<chunksize=(1, 120, 50, 10, 8), meta=np.ndarray>
    Tz         (node, time, k_l, j, i) float32 dask.array<chunksize=(1, 120, 1, 10, 8), meta=np.ndarray>
    GM_CHI     (node, time, k, j, i) float32 dask.array<chunksize=(1, 120, 49, 10, 8), meta=np.ndarray>
Attributes:
    Conventions:  CF-1.6
    title:        netCDF wrapper of MITgcm MDS binary data
    source:       MITgcm
    history:      Created by calling `open_mdsdataset(grid_dir='/glade/work/d...
(natre,) = dask.optimize(natre)

natre_mean = natre.cf.mean(["X", "Y"]).load()

/glade/u/home/dcherian/miniconda3/envs/eddydiff/lib/python3.10/site-packages/dask_jobqueue/core.py:255: FutureWarning: job_extra has been renamed to job_extra_directives. You are still using it (even if only set to []; please also check config files). If you did not set job_extra_directives yet, job_extra will be respected for now, but it will be removed in a future release. If you already set job_extra_directives, job_extra is ignored and you can remove it.
  warnings.warn(warn, FutureWarning)
/glade/u/home/dcherian/miniconda3/envs/eddydiff/lib/python3.10/site-packages/dask_jobqueue/core.py:274: FutureWarning: env_extra has been renamed to job_script_prologue. You are still using it (even if only set to []; please also check config files). If you did not set job_script_prologue yet, env_extra will be respected for now, but it will be removed in a future release. If you already set job_script_prologue, env_extra is ignored and you can remove it.
  warnings.warn(warn, FutureWarning)

natre_mean.to_netcdf("../datasets/ecco-chi.nc")

natre_mean["GM_CHIX"].mean("time").cf.plot(hue="node")

[<matplotlib.lines.Line2D at 0x2ac11e83d180>,
 <matplotlib.lines.Line2D at 0x2ac11e83d030>]
_images/5a6bb9edf2cb49b016f07ac0969a327fe4b16102d7dcc5f4e56cd400789be731.png 
natre_mean["GM_Kux"].mean("time").cf.plot(hue="node")
natre_mean["GM_Kvy"].mean("time").cf.plot(hue="node")

[<matplotlib.lines.Line2D at 0x2abef70a2380>,
 <matplotlib.lines.Line2D at 0x2abef70a38b0>]
_images/536651fcf4b0fe54c3b4b0b5f8053f1885f3b8908af5ef8f5dcd8cb20e487644.png 
natre_mean.GM_Kvy.mean("time").cf.plot(hue="node")

[<matplotlib.lines.Line2D at 0x2abe4ea55c90>,
 <matplotlib.lines.Line2D at 0x2abe4ea57160>]
_images/fee7881ae3c7d3d023ce5927ae83c486308528841f3839cc4b814f07e4a2ce55.png

Dev#

dsfull, grid = open_directory("/glade/work/dcherian/mitgcm/ECCOV4/release4/run/diags/")
dsfull

ds = dsfull.sel(face=2).load()

/glade/u/home/dcherian/miniconda3/envs/eddydiff/lib/python3.10/site-packages/dask_jobqueue/core.py:255: FutureWarning: job_extra has been renamed to job_extra_directives. You are still using it (even if only set to []; please also check config files). If you did not set job_extra_directives yet, job_extra will be respected for now, but it will be removed in a future release. If you already set job_extra_directives, job_extra is ignored and you can remove it.
  warnings.warn(warn, FutureWarning)
/glade/u/home/dcherian/miniconda3/envs/eddydiff/lib/python3.10/site-packages/dask_jobqueue/core.py:274: FutureWarning: env_extra has been renamed to job_script_prologue. You are still using it (even if only set to []; please also check config files). If you did not set job_script_prologue yet, env_extra will be respected for now, but it will be removed in a future release. If you already set job_script_prologue, env_extra is ignored and you can remove it.
  warnings.warn(warn, FutureWarning)

NATRE region#

natre = ds.isel(i=slice(10, 18), j=slice(10, 20), i_g=slice(10, 18), j_g=slice(10, 20))

natre.load(scheduler=client)

<xarray.Dataset>
Dimensions:    (i: 8, i_g: 8, j: 10, j_g: 10, k: 50, k_u: 50, k_l: 50,
                k_p1: 51, time: 24)
Coordinates: (12/45)
  * i          (i) int64 10 11 12 13 14 15 16 17
  * i_g        (i_g) int64 10 11 12 13 14 15 16 17
  * j          (j) int64 10 11 12 13 14 15 16 17 18 19
  * j_g        (j_g) int64 10 11 12 13 14 15 16 17 18 19
  * k          (k) int64 0 1 2 3 4 5 6 7 8 9 ... 40 41 42 43 44 45 46 47 48 49
  * k_u        (k_u) int64 0 1 2 3 4 5 6 7 8 9 ... 40 41 42 43 44 45 46 47 48 49
    ...         ...
    maskCtrlS  (k, j_g, i) bool True True True True ... False False False False
    rhoRef     (k) >f4 1.024e+03 1.024e+03 1.024e+03 ... 1.052e+03 1.054e+03
    maskCtrlW  (k, j, i_g) bool True True True True ... False False False False
    iter       (time) int64 1 2 3 4 5 6 7 8 9 10 ... 16 17 18 19 20 21 22 23 24
  * time       (time) timedelta64[ns] 00:00:01 00:00:02 ... 00:00:23 00:00:24
    drC_kl     (k_l) >f4 10.0 10.0 10.0 10.0 10.0 ... 399.0 422.0 445.0 228.2
Data variables: (12/19)
    GM_CHIZ    (time, k_l, j, i) float32 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0
    GM_ubT     (time, k, j, i_g) float32 1.363e+06 1.334e+06 ... 0.0 0.0
    GM_CHIX    (time, k, j, i_g) float32 1.571e-09 3.038e-09 ... 0.0 0.0
    GM_Kvz     (time, k, j_g, i) float32 3.523 3.411 3.31 3.206 ... 0.0 0.0 0.0
    GM_Kwz     (time, k_l, j, i) float32 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0
    THETA      (time, k, j, i) float32 23.1 22.98 22.81 22.6 ... 0.0 0.0 0.0 0.0
    ...         ...
    Tx         (time, k, j, i_g) float32 -8.891e-07 -1.213e-06 ... 0.0 0.0
    Ty         (time, k, j_g, i) float32 -2.024e-06 -1.976e-06 ... 0.0 0.0
    Tz         (time, k_l, j, i) float32 nan nan nan nan nan ... 0.0 0.0 0.0 0.0
    chix_      (time, k, j, i_g) float32 1.433e-09 2.585e-09 ... 0.0 0.0
    chiy_      (time, k, j_g, i) float32 7.264e-09 6.708e-09 ... 0.0 0.0
    GM_CHI     (time, k, j, i) float32 1.467e-08 2.126e-08 2.9e-08 ... nan nan
Attributes:
    Conventions:  CF-1.6
    title:        netCDF wrapper of MITgcm MDS binary data
    source:       MITgcm
    history:      Created by calling `open_mdsdataset(grid_dir='/glade/work/d...
natre.GM_CHI.mean(["time", "j", "i"]).plot(y="Z")

[<matplotlib.lines.Line2D at 0x2b4fc22cb430>]
_images/f62bdb531d483187e4c5c0be37188ac2da7f383b4dceebf31fea607bd340ae81.png

Development#

Checking CHIX, CHIXO, CHIY, CHIYO#

ds["chix_"] = ds.GM_Kux * ds.Tx**2
ds["chixo_"] = 2 * ds.GM_Kuz * ds.Tx * grid.interp(ds.Tz, axis=["X", "Z"])
ds["chiy_"] = ds.GM_Kvy * ds.Ty**2
ds["chiyo_"] = 2 * ds.GM_Kvz * ds.Ty * grid.interp(ds.Tz, axis=["Y", "Z"])

ds[["GM_Kvz"]].isel(time=5, j_g=20, i=20).to_array().plot(hue="variable")
ds[["GM_Kuz"]].isel(time=5, j=20, i_g=20).to_array().plot(hue="variable")

[<matplotlib.lines.Line2D at 0x2ab5eb6bc550>]
_images/4390a5c4b063bd90c2faa0dd09ec52d95a87584474d96f7ee2fbcbb2d389bef3.png 
ds[["chix_", "GM_CHIX"]].isel(time=5, j=20, i_g=20).to_array().plot(
    hue="variable", yscale="log"
)

[<matplotlib.lines.Line2D at 0x2ab5f2fcf730>,
 <matplotlib.lines.Line2D at 0x2ab5f2fcf700>]
_images/29c3302323ffc820cac14b5ff4a3fa1aacabba4bd49d45cb2887bb31f4dd65da.png 
ds[["chixo_", "GM_CHIXO"]].isel(time=5, j=20, i_g=20).to_array().plot(hue="variable")

[<matplotlib.lines.Line2D at 0x2ab5eee5f400>,
 <matplotlib.lines.Line2D at 0x2ab5eee5e1a0>]
_images/faebae1ed55323ad9517ffcbd7239fc1f9a9b9ab4938f35ecbec79caee80dd50.png 
ds[["chiy_", "GM_CHIY"]].isel(time=5, j_g=20, i=20).to_array().plot(
    hue="variable", yscale="log"
)

[<matplotlib.lines.Line2D at 0x2ab60a3ef8e0>,
 <matplotlib.lines.Line2D at 0x2ab60a3ee0e0>]
_images/5d0c3308367df0c08a0c09fedb86fd96aa47d4ff3317485789f7be35e51437bd.png 
ds[["chiyo_", "GM_CHIYO"]].isel(time=5, j_g=20, i=20).to_array().plot(
    hue="variable",
)

[<matplotlib.lines.Line2D at 0x2ab5f33cbac0>,
 <matplotlib.lines.Line2D at 0x2ab5f33cba90>]
_images/7f5b6d891427f7db82dc8a9be3486ac72c0399b63e2fa5dc195d352d4415c8d1.png

Checking CHIZ#

subset = ds.isel(time=4, i_g=4, i=4, j=39, j_g=39)
subset["GM_Tz"] = subset["GM_KwzTz"] / subset.GM_Kwz / subset.rA
# subset["THETA"] = ds.THETA.isel(time=5,i=4, j=39)
subset = subset.load()

subset.drC.isel(k_p1=slice(1, None)) still seems right

subset["Tz_2"] = (
    -grid.diff(subset.THETA, "Z")
    / subset.drC.isel(k_p1=slice(1, None)).rename({"k_p1": "k_l"}).variable
)

# (subset["GM_KwzTz"] * subset["GM_KwzTz"] / subset.GM_Kwz / subset.rA**2).plot()
np.sqrt(subset.GM_CHIZ / subset.GM_Kwz).plot(marker=".")
np.abs(subset.Tz).plot()
np.abs(subset.Tz_2).plot()

[<matplotlib.lines.Line2D at 0x2ab52a7568c0>]
_images/814d01d4cbfdbd902c75b9cf09977a07e6835f009e01061f1ce13ddb61b3fc8c.png 
(subset["GM_KwzTz"]).plot(marker=".")
(-subset.GM_Kwz * subset.Tz * subset.rA).plot()
(-subset.GM_Kwz * subset.Tz_2 * subset.rA).plot()

[<matplotlib.lines.Line2D at 0x2ab5f4135660>]
_images/dab564147375c8c7331b3759303ebc6fda2fce22390710b86436bf524119c9eb.png

Is Kvz == Kwy?#

Seems close, and errors from averaging probably?

grid.interp(ds.GM_Kvz, axis="Z").isel(j_g=20, i=20, time=1).plot()
grid.interp(ds.GM_Kwy, axis="Y").isel(j_g=20, i=20, time=1).plot()

[<matplotlib.lines.Line2D at 0x2b4fc513cc70>]
_images/07e0d8d7acb1e5b17a882c4f7a58c146ce983b1f82a27234f2c0f17c4ef25633.png

File reading experiments#

ecco.read_llc_to_tiles(
    dirname + "GM_CHIZ_INST/", "GM_CHIZ.0000000001.data", nk=50, nl=1, use_xmitgcm=True
)

read_llc_to_tiles: full_filename:  /glade/work/dcherian/mitgcm/ECCOV4/release4/run/diags/GM_CHIZ_INST//GM_CHIZ.0000000001.data
Array Chunk
Bytes 20.08 MiB 31.64 kiB
Shape (1, 50, 13, 90, 90) (1, 1, 1, 90, 90)
Dask graph 650 chunks in 1 graph layer
Data type >f4 numpy.ndarray
50 1 90 90 13 
import xmitgcm

xmitgcm.open_mdsdataset(dirname + "GM_CHIZ_INST/", dirname + "../")

<xarray.Dataset>
Dimensions:    (XC: 90, YC: 1170, XG: 90, YG: 1170, Z: 50, Zp1: 51, Zu: 50,
                Zl: 50, time: 8)
Coordinates: (12/33)
  * XC         (XC) >f4 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0 0.0
  * YC         (YC) >f4 0.0 0.0 0.0 0.0 0.0 ... 9.482 -57.27 67.47 9.482 -57.27
  * XG         (XG) >f4 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0 0.0
  * YG         (YG) >f4 0.0 0.0 0.0 0.0 0.0 ... 9.97 -57.01 67.56 9.97 -57.01
  * Z          (Z) >f4 -5.0 -15.0 -25.0 ... -5.039e+03 -5.461e+03 -5.906e+03
  * Zp1        (Zp1) >f4 0.0 -10.0 -20.0 ... -5.244e+03 -5.678e+03 -6.134e+03
    ...         ...
    maskCtrlC  (Z, YC, XC) bool dask.array<chunksize=(50, 1170, 90), meta=np.ndarray>
    maskCtrlS  (Z, YG, XC) bool dask.array<chunksize=(50, 1170, 90), meta=np.ndarray>
    rhoRef     (Z) >f4 dask.array<chunksize=(50,), meta=np.ndarray>
    maskCtrlW  (Z, YC, XG) bool dask.array<chunksize=(50, 1170, 90), meta=np.ndarray>
    iter       (time) int64 dask.array<chunksize=(1,), meta=np.ndarray>
  * time       (time) timedelta64[ns] 00:00:01 00:00:02 ... 00:00:07 00:00:08
Data variables:
    GM_CHIZ    (time, Z, YC, XC) float32 dask.array<chunksize=(1, 50, 1170, 90), meta=np.ndarray>
Attributes:
    Conventions:  CF-1.6
    title:        netCDF wrapper of MITgcm MDS binary data
    source:       MITgcm
    history:      Created by calling `open_mdsdataset(grid_dir='/glade/work/d...
chiz = ecco.llc_tiles_to_xda(
    ecco.read_llc_to_tiles(
        dirname + "GM_CHIZ_INST/",
        "GM_CHIZ.0000000001.data",
        nk=50,
        nl=1,
        use_xmitgcm=True,
    ),
    var_type="c",
    dim4="depth",
    dim5="time",
)
chiz

read_llc_to_tiles: full_filename:  /glade/work/dcherian/mitgcm/ECCOV4/release4/run/diags/GM_CHIZ_INST//GM_CHIZ.0000000001.data

<xarray.DataArray 'llc-2d8997a5c3544544991ec6632d79eee9' (time: 1, k: 50,
                                                          tile: 13, j: 90, i: 90)>
dask.array<llc, shape=(1, 50, 13, 90, 90), dtype=>f4, chunksize=(1, 1, 1, 90, 90), chunktype=numpy.ndarray>
Coordinates:
  * k        (k) int64 0 1 2 3 4 5 6 7 8 9 10 ... 40 41 42 43 44 45 46 47 48 49
  * time     (time) int64 0
  * tile     (tile) int64 0 1 2 3 4 5 6 7 8 9 10 11 12
  * j        (j) int64 0 1 2 3 4 5 6 7 8 9 10 ... 80 81 82 83 84 85 86 87 88 89
  * i        (i) int64 0 1 2 3 4 5 6 7 8 9 10 ... 80 81 82 83 84 85 86 87 88 89
chix.sel(tile=5).plot()

<matplotlib.collections.QuadMesh at 0x2b9eaf65a830>
_images/3f047c710324b9986afabee88c3cf101b3d4263621be202047fb2e967dae7d6f.png

Tracking down gradient calculation in gmredi_ytransport.F#

Here, I saved the dTdz in gmredi_ytransport.F as GM_CHIYO

dTdz_v replicates the calculation.

This helped tracked down a bug in ds.drC_kl (we want drC.isel(kp1=slice(-1)) instead of drC.isel(kp1=slice(1, None))

ΔT_z = -1 * grid.diff(ds.THETA, "Z")
ΔT_z_y = grid.interp(ΔT_z, "Y")
dTdz_v = grid.interp(ΔT_z_y / ds.drC_kl, "Z")

(ds.GM_CHIYO).isel(time=5, j_g=21, i=20).plot()
(dTdz_v).isel(time=5, j_g=21, i=20).plot()
(grid.interp(ds.Tz, axis=["Z", "Y"])).isel(time=5, j_g=21, i=20).plot()

[<matplotlib.lines.Line2D at 0x2ab5ea93fd00>]
_images/5432cca27df52cf3a69e3220c745c08dc1449528bd2b20d5024bb254a07009da.png

Test ECCO Kredi application#

In input_init/README(link) we see

    total_diffkr_r009bit11.bin           vert. diff. of release 1 (this field plus xx is the total)
    total_kapgm_r009bit11.bin            Kappa GM of release 1 (this field plus xx is the total)
    total_kapredi_r009bit11.bin          Kappa Redi of release 1 (this field plus xx is the total)
    xx_*.*                               control adjustments

I suspect by commenting out the file name, I set a constant and the control adjustments were applied

And it turns out the control adjustments need to be rescaled so that makes it a pain to check.

It does approximately match so that seems to be what’s going on.

path = "/glade/u/home/dcherian/work/mitgcm/ECCOV4/release4/input_init/"

import xarray as xr

natre_mean = xr.open_dataset("../datasets/ecco-chi.nc").isel(time=1)

import numpy as np
import xmitgcm

kapredi_adj = xmitgcm.open_mdsdataset(
    path,
    grid_dir=path + "../run/",
    prefix="xx_kapredi",
    geometry="llc",
    extra_variables=dict(xx_kapredi=dict(dims=["k", "j", "i"], attrs={})),
).load()
kapredi_adj

<xarray.Dataset>
Dimensions:     (i: 90, i_g: 90, j: 90, j_g: 90, k: 50, k_u: 50, k_l: 50,
                 k_p1: 51, face: 13, time: 1)
Coordinates: (12/44)
  * i           (i) int64 0 1 2 3 4 5 6 7 8 9 ... 80 81 82 83 84 85 86 87 88 89
  * i_g         (i_g) int64 0 1 2 3 4 5 6 7 8 9 ... 81 82 83 84 85 86 87 88 89
  * j           (j) int64 0 1 2 3 4 5 6 7 8 9 ... 80 81 82 83 84 85 86 87 88 89
  * j_g         (j_g) int64 0 1 2 3 4 5 6 7 8 9 ... 81 82 83 84 85 86 87 88 89
  * k           (k) int64 0 1 2 3 4 5 6 7 8 9 ... 40 41 42 43 44 45 46 47 48 49
  * k_u         (k_u) int64 0 1 2 3 4 5 6 7 8 9 ... 41 42 43 44 45 46 47 48 49
    ...          ...
    maskW       (k, face, j, i_g) bool False False False ... False False False
    maskS       (k, face, j_g, i) bool False False False ... False False False
    maskCtrlC   (k, face, j, i) bool False False False ... False False False
    maskCtrlS   (k, face, j_g, i) bool False False False ... False False False
    rhoRef      (k) >f4 1.024e+03 1.024e+03 1.024e+03 ... 1.052e+03 1.054e+03
    maskCtrlW   (k, face, j, i_g) bool False False False ... False False False
Data variables:
    xx_kapredi  (time, k, face, j, i) >f4 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0
Attributes:
    Conventions:  CF-1.6
    title:        netCDF wrapper of MITgcm MDS binary data
    source:       MITgcm
    history:      Created by calling `open_mdsdataset(grid_dir='/glade/u/home...
natre_mean.GM_Kux.load()

<xarray.DataArray 'GM_Kux' (node: 2, k: 50)>
array([[1653.3713 , 1756.2273 , 1774.7654 , 1698.0707 , 1606.0331 ,
        1502.6423 , 1411.8212 , 1353.8236 , 1326.3242 , 1335.0173 ,
        1382.8555 , 1460.9812 , 1550.7936 , 1637.4607 , 1704.886  ,
        1736.8871 , 1719.9307 , 1662.3121 , 1585.8231 , 1500.4946 ,
        1408.286  , 1299.3516 , 1193.2825 , 1113.6604 , 1086.0511 ,
        1114.4277 , 1220.7032 , 1396.9984 , 1558.5225 , 1537.3969 ,
        1229.34   ,  750.17395,  366.39853,  213.548  ,  228.79117,
         287.18646,  343.04477,  408.48822,  487.57626,  579.28723,
         710.1041 ,  898.3161 , 1090.6859 , 1212.5089 , 1229.1472 ,
        1170.0338 , 1088.2152 , 1028.699  ,  999.34344, 1000.     ],
       [1122.2953 , 1138.5017 , 1138.5872 , 1127.2064 , 1121.2599 ,
        1119.6086 , 1118.9813 , 1114.1803 , 1100.417  , 1080.9762 ,
        1060.82   , 1044.8038 , 1036.1597 , 1035.7036 , 1040.701  ,
        1040.6842 , 1025.8724 ,  993.8615 ,  949.4629 ,  908.0531 ,
         887.2899 ,  903.4036 ,  956.21436, 1028.5127 , 1086.1965 ,
        1094.7805 , 1055.2051 ,  999.22345,  958.0874 ,  941.8081 ,
         944.84357,  962.67786,  988.8215 , 1001.7471 ,  980.7318 ,
         939.7799 ,  920.5631 ,  944.925  ,  995.82355, 1041.636  ,
        1061.2731 , 1056.8977 , 1041.72   , 1026.6104 , 1016.8036 ,
        1011.05994, 1007.35706, 1004.25977, 1000.74286, 1000.     ]],
      dtype=float32)
Coordinates:
  * k        (k) int64 0 1 2 3 4 5 6 7 8 9 10 ... 40 41 42 43 44 45 46 47 48 49
    face     int64 2
    time     timedelta64[ns] 00:23:48
    Z        (k) float32 -5.0 -15.0 -25.0 ... -5.039e+03 -5.461e+03 -5.906e+03
    drF      (k) float32 10.0 10.0 10.0 10.0 10.0 ... 387.5 410.5 433.5 456.5
    PHrefC   (k) float32 49.05 147.1 245.2 ... 4.944e+04 5.357e+04 5.794e+04
    rhoRef   (k) float32 1.024e+03 1.024e+03 1.024e+03 ... 1.052e+03 1.054e+03
    iter     float64 1.428e+03
  * node     (node) object 'adj' 'const'
(kapredi_adj).sel(face=2).isel(time=0, i=slice(10, 18), j=slice(10, 20)).mean(
    ["i", "j"]
).xx_kapredi.plot(color="k")
(natre_mean.GM_Kux / 2000 - 0.5).sel(node="const").plot.line(hue="node")

[<matplotlib.lines.Line2D at 0x2ba4c4f1bd60>]
_images/2902c3acec4b8700edd0e4e2fa6b1b858e3bf2fa2e780b876fe6a19a01ae16ee.png