LES tendencies#

  1. This has to be run on casper

  2. I would add a 5 day buffer on les_time_length

  3. I’ve removed the coriolis contribution, so the LES needs f≠0

  4. I have not checked the 12H offset in time

/glade/p/cgd/oce/people/dwhitt/TPOS/tpos_LES_runs_setup_scripts/tpos-DIABLO/diablo_2.0/tpos_config_files

import glob

import cf_xarray
import dask
import dcpy
import h5py
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import xarray as xr
import xgcm

import pump
import pump.les
import pump.model
from pump.les import interleave_time, interpolate, write_to_txt

dask.config.set(**{"array.slicing.split_large_chunks": False})  # silence some warnings
<dask.config.set at 0x2b963be83ee0>

Sets up a dask cluster; only run once

import ncar_jobqueue

cluster = ncar_jobqueue.NCARCluster(account="ncgd0011")
cluster.scale(12)
import distributed

client = distributed.Client(cluster)
client

Client

Cluster

  • Workers: 0
  • Cores: 0
  • Memory: 0 B

Things you could change; if changing, rerun all cells below (Use the Run menu above)

sst = xr.open_mfdataset(
    f"{gcmdir}/*buoy*.nc",
    chunks={"latitude": 120, "longitude": 500, "depth": 1},
    combine="by_coords",
    parallel=True,
).theta.isel(depth=0)
f, axx = plt.subplots(10, 2, sharey=True, constrained_layout=True)
for year, ax in zip(np.unique(sst.time.dt.year), axx.flat):
    (
        sst.theta.isel(depth=0)
        .sel(latitude=slice(-2, 5))
        .sel(time=slice(f"{year}-Jun", f"{year}-Dec"))
        .sel(longitude=-110, method="nearest")
        .plot(
            x="time",
            ax=ax,
            vmin=22,
            vmax=28,
            cmap=mpl.cm.RdYlBu_r,
            add_colorbar=False,
            add_labels=False,
        )
    )

f.set_size_inches((12, 18))
_images/7f11aa6ff952bb4b11637928be0a2bc841ba4780903ab71145dc12c079af521e.png

Choose time period + plot ocean state to see what we get#

gcmdir = "/glade/campaign/cgd/oce/people/bachman/TPOS_1_20_20_year_output/"  # MITgcm output directory
outdir = "/glade/work/dcherian/pump/les_forcing/mitgcm_test/"  # output directory for txt files

whichz = "216"  # or "288"  # which z grid h5 file to use

leslon = -110  # it will find the nearest
leslat = 0  # will find the nearest

# start date for les; noon is a good time (it is before sunrise)
sim_time = pd.Timestamp("2015-09-20 12:00:00")

# ADD a 5 day buffer here (:] all kinds of bugs at the beginning and end)
les_time_length = 30  # (days); length of time for forcing/pushing files

# don't change anything here
output_start_time = pd.Timestamp("1999-01-01")  # don't change
firstfilenum = (
    (
        sim_time - pd.Timedelta("1D") - output_start_time
    )  # add one day offset to avoid bugs
    .to_numpy()
    .astype("timedelta64[D]")
    .astype("int")
)
lastfilenum = firstfilenum + les_time_length


def gen_file_list(suffix):
    files = [
        f"{gcmdir}/File_{num}_{suffix}.nc"
        for num in range(firstfilenum - 1, lastfilenum + 1)
    ]

    return files


# station = pump.model.model.read_stations_20(
#    dayglobstr=folder, dirname=gcmdir, globstr="*_-140*.nc"
# )

state = xr.open_mfdataset(
    gen_file_list(suffix="buoy"),
    chunks={"latitude": 120, "longitude": 500},
    combine="by_coords",
    parallel=True,
).rename({"latitude": "YC", "longitude": "XC", "depth": "RC"})
state["YC"].attrs["units"] = "degrees_north"


q, r = divmod(firstfilenum, 250)
folder = f"{250*q + 1}-{250*(q+1)}"
path = f"{gcmdir}/STATION_DATA/Day_{folder}/"

station = (
    xr.open_mfdataset(
        sorted(glob.glob(path + f"/*_{leslon}.000.nc")),
        parallel=True,
        decode_times=False,
    )
    .sel(depth=slice(-150))
    .squeeze()
)
station.time.attrs["units"] = f"seconds since {output_start_time}"
station = xr.decode_cf(station)
station["dens"] = pump.mdjwf.dens(station.salt, station.theta, np.array([0]))
station["mld"] = pump.calc.get_mld(station.dens).compute()
station

choose_lats = [-0.25, 2]
subset = station.sel(time=slice(sim_time, sim_time + pd.Timedelta("30D"))).sel(
    latitude=choose_lats, method="nearest"
)
subset["sst"] = state.theta.sel(XC=leslon, method="nearest", RC=0).reindex(
    time=subset.time, method="pad"
)
subset = pump.calc.calc_reduced_shear(subset)
subset["eucmax"] = pump.calc.get_euc_max(subset.isel(latitude=0).u)
subset["dcl_base"] = pump.calc.get_dcl_base_Ri(subset)
subset["Jq"] = 1035 * 3995 * subset.DFrI_TH / 3.09e7

f, ax = pump.plot.plot_jq_sst(
    subset,
    full=subset,
    lon=leslon,
    lat=choose_lats,
    periods=slice(subset.time[0].values, subset.time[-1].values),
)
ax[2, 0].set_ylim([-150, 0])
calc uz
calc vz
calc S2
calc N2
calc shred2
Calc Ri
(-150.0, 0.0)
_images/5785a5c305a3e47b50272a39bf8deab0e1b7e4e58a819341cb2f6d1e442f271b.png
state.theta.sel(XC=leslon, method="nearest", RC=0).plot(
    x="time", robust=True, cmap=mpl.cm.RdYlBu_r
)
dcpy.plots.liney(choose_lats, zorder=1)

fg = subset.DFrI_TH.sortby("latitude", ascending=False).plot(
    row="latitude", robust=True, cbar_kwargs={"orientation": "horizontal"}
)
[
    (-1 * subset.KPPhbl).sel(loc).plot(color="r", ax=ax, zorder=12, _labels=False)
    for loc, ax in zip(fg.name_dicts.squeeze(), fg.axes.squeeze())
]
[
    (subset.mld).sel(loc).plot(color="C1", ax=ax, zorder=13, _labels=False)
    for loc, ax in zip(fg.name_dicts.squeeze(), fg.axes.squeeze())
]

fg.fig.set_size_inches((10, 5))
_images/6509b0d1a0f6dcf56aeb16c89f27b4ba5e2e11c2aa3dafabb7551a33db1394cc.png _images/51ad1534a8d2437ab8cb5273f91161bdfce45bc2704013e412b51afb7f2acd3f.png

LES config#

LES things that you need not change; Grid is chosen by whichz in the above cell

# Should not need to change anything below this
griddir = "/glade/p/cgd/oce/people/dwhitt/TPOS/tpos_LES_runs_setup_scripts/tpos-DIABLO/diablo_2.0/pre_process/"
renamer = {
    "u": "U",
    "v": "W",
    "temp": "T",
    # "salt": "S"
}
offset = xr.Dataset({"U": 0.4, "T": 25, "S": 35.25, "W": 0})
t0 = (sim_time).to_numpy()

grid = h5py.File(f"{griddir}/gridtpos{whichz}.h5", mode="r")
zles = grid["grids"]["y"][:]
grid.close()
zles = 0.5 * (zles[1:-1] + zles[2:])
newz = xr.DataArray(
    zles - zles.max() - 0.25, dims="z", name="z"  # this is what dan does
)
newz  # prints out the grid
<xarray.DataArray 'z' (z: 216)>
array([-107.75, -107.25, -106.75, -106.25, -105.75, -105.25, -104.75,
       -104.25, -103.75, -103.25, -102.75, -102.25, -101.75, -101.25,
       -100.75, -100.25,  -99.75,  -99.25,  -98.75,  -98.25,  -97.75,
        -97.25,  -96.75,  -96.25,  -95.75,  -95.25,  -94.75,  -94.25,
        -93.75,  -93.25,  -92.75,  -92.25,  -91.75,  -91.25,  -90.75,
        -90.25,  -89.75,  -89.25,  -88.75,  -88.25,  -87.75,  -87.25,
        -86.75,  -86.25,  -85.75,  -85.25,  -84.75,  -84.25,  -83.75,
        -83.25,  -82.75,  -82.25,  -81.75,  -81.25,  -80.75,  -80.25,
        -79.75,  -79.25,  -78.75,  -78.25,  -77.75,  -77.25,  -76.75,
        -76.25,  -75.75,  -75.25,  -74.75,  -74.25,  -73.75,  -73.25,
        -72.75,  -72.25,  -71.75,  -71.25,  -70.75,  -70.25,  -69.75,
        -69.25,  -68.75,  -68.25,  -67.75,  -67.25,  -66.75,  -66.25,
        -65.75,  -65.25,  -64.75,  -64.25,  -63.75,  -63.25,  -62.75,
        -62.25,  -61.75,  -61.25,  -60.75,  -60.25,  -59.75,  -59.25,
        -58.75,  -58.25,  -57.75,  -57.25,  -56.75,  -56.25,  -55.75,
        -55.25,  -54.75,  -54.25,  -53.75,  -53.25,  -52.75,  -52.25,
        -51.75,  -51.25,  -50.75,  -50.25,  -49.75,  -49.25,  -48.75,
        -48.25,  -47.75,  -47.25,  -46.75,  -46.25,  -45.75,  -45.25,
        -44.75,  -44.25,  -43.75,  -43.25,  -42.75,  -42.25,  -41.75,
        -41.25,  -40.75,  -40.25,  -39.75,  -39.25,  -38.75,  -38.25,
        -37.75,  -37.25,  -36.75,  -36.25,  -35.75,  -35.25,  -34.75,
        -34.25,  -33.75,  -33.25,  -32.75,  -32.25,  -31.75,  -31.25,
        -30.75,  -30.25,  -29.75,  -29.25,  -28.75,  -28.25,  -27.75,
        -27.25,  -26.75,  -26.25,  -25.75,  -25.25,  -24.75,  -24.25,
        -23.75,  -23.25,  -22.75,  -22.25,  -21.75,  -21.25,  -20.75,
        -20.25,  -19.75,  -19.25,  -18.75,  -18.25,  -17.75,  -17.25,
        -16.75,  -16.25,  -15.75,  -15.25,  -14.75,  -14.25,  -13.75,
        -13.25,  -12.75,  -12.25,  -11.75,  -11.25,  -10.75,  -10.25,
         -9.75,   -9.25,   -8.75,   -8.25,   -7.75,   -7.25,   -6.75,
         -6.25,   -5.75,   -5.25,   -4.75,   -4.25,   -3.75,   -3.25,
         -2.75,   -2.25,   -1.75,   -1.25,   -0.75,   -0.25])
Dimensions without coordinates: z

MITgcm#

The daily avg output was subset in MATLAB syntax : 41:end-40, 41:end-40, 1:nz

Read all ocean model output#

coords = pump.model.read_mitgcm_coords(gcmdir)
grid = xgcm.Grid(
    coords, periodic=False, boundary={"X": "extend", "Y": "extend", "Z": "extend"}
)

# ssh data
ssh = (
    xr.open_mfdataset(
        gen_file_list(suffix="etan"),
        chunks={"latitude": 120, "longitude": 500},
        combine="by_coords",
        parallel=True,
    )
    .rename({"latitude": "YC", "longitude": "XC"})
    .update(coords.coords)
)

# momentum budget terms
budgets = xr.open_mfdataset(
    gen_file_list(suffix="ub") + gen_file_list(suffix="vb"),
    chunks={"latitude": 120, "longitude": 500},
    combine="by_coords",
    parallel=True,
)
budgets = pump.model.rename_mitgcm_budget_terms(budgets, coords).update(coords.coords)

# surface data
surf = (
    xr.open_mfdataset(
        gen_file_list(suffix="surf"),
        chunks={"latitude": 120, "longitude": 500},
        combine="by_coords",
        parallel=True,
    )
    .rename({"latitude": "YC", "longitude": "XC"})
    .update(coords.coords)
)
surf["oceQsw"] = surf.oceQsw.fillna(0)
surf["taux"] = 1035 * 2.5 * budgets["Um_Ext"].isel(RC=0)
surf["tauy"] = 1035 * 2.5 * budgets["Vm_Ext"].isel(RC=0)

# heat budget terms
hb = xr.open_mfdataset(
    gen_file_list(suffix="hb"),
    chunks={"latitude": 120, "longitude": 500},
    combine="by_coords",
    parallel=True,
)
hb_ren = pump.model.rename_mitgcm_budget_terms(hb, coords).update(coords.coords)

# state variables
state = xr.open_mfdataset(
    gen_file_list(suffix="buoy"),
    chunks={"latitude": 120, "longitude": 500},
    combine="by_coords",
    parallel=True,
).rename({"latitude": "YC", "longitude": "XC", "depth": "RC"})
state["u"] = state.u.rename({"XC": "XG"})
state["v"] = state.v.rename({"YC": "YG"})
state = state.update(coords.coords)


# grid metrics
metrics = (
    pump.model.read_metrics(gcmdir)
    .compute()
    .isel(longitude=slice(40, -40), latitude=slice(40, -40), depth=slice(136))
    .chunk({"latitude": 120, "longitude": 500})
    .drop(["latitude", "longitude", "depth"])
    .rename({"depth_left": "RF", "depth": "RC", "latitude": "YC", "longitude": "XC"})
    .set_coords("RF")
    .isel(RF=slice(136))
)
metrics["rAw"] = metrics.rAw.rename({"XC": "XG"})
metrics["rAs"] = metrics.rAs.rename({"YC": "YG"})
metrics["hFacW"] = metrics.hFacW.rename({"XC": "XG"})
metrics["hFacS"] = metrics.hFacS.rename({"YC": "YG"})
metrics["DXG"] = metrics.DXG.rename({"YC": "YG"})
metrics["DYG"] = metrics.DYG.rename({"XC": "XG"})
metrics = metrics.update(coords.coords)

# penetrative shortwave radiation profile
swprofile = (
    xr.DataArray(
        pump.model.mitgcm_sw_prof(metrics.RF[:-1].values)
        - pump.model.mitgcm_sw_prof(metrics.RF[1:].values),
        dims=["RC"],
        attrs={"long_name": "SW radiation deposited in cell"},
    )
    .assign_coords(RC=coords.coords["RC"][:-1])
    .reindex(RC=coords.coords["RC"])
)

Offset in time (I have not checked this).

ssh["time"] = ssh.time + pd.Timedelta("12H")
budgets["time"] = budgets.time + pd.Timedelta("12H")
state["time"] = state.time + pd.Timedelta("12h")
hb_ren["time"] = hb_ren.time + pd.Timedelta("12h")
surf["time"] = surf.time + pd.Timedelta("12h")

Make sure that everything is on the same grid; just a check.

xr.align(budgets, ssh, state, hb_ren, metrics, surf, join="exact");

Estimate flux divergences#

# constants taken from diagnostics.log (I think)
global_area = 2.196468634481708e13
rhoConst = 1035
Cp = 3994
f0 = 2 * (2 * np.pi / 86400) * np.sin(leslat * np.pi / 180)

# This is WRONG! Because we have subset the output this won't recover the MITgcm correction exactly
# seems to be particularly bad at t=0
surf_mass = (hb_ren.WTHMASS * metrics.RAC).isel(RC=0)
TsurfCor = surf_mass.sum(["XC", "YC"]) / global_area


def expand_depth(ds):
    if "RC" not in ds:
        ds["RC"] = coords.RC.isel(RC=0).data
    return ds.expand_dims("RC").reindex(RC=coords.RC, fill_value=0)


budgets["surf_corr_tend"] = -1 * (
    (TsurfCor - hb_ren.WTHMASS.isel(RC=0))
    / metrics.dRF.isel(RC=0)
    / metrics.hFacC.isel(RC=0)
).pipe(expand_depth)
budgets["TOTTTEND"] = hb_ren.TOTTTEND


budgets["ADV_TH"] = (
    -1
    * (
        grid.diff(hb_ren.ADVx_TH, "X")
        + grid.diff(hb_ren.ADVy_TH, "Y")
        - grid.diff(
            hb_ren.ADVr_TH.fillna(0), "Z"
        )  # negative because level 0 is surface
    )
    / metrics.cellvol
)

budgets["MIX_TH"] = (
    grid.diff(hb_ren.DFrI_TH.fillna(0) + hb_ren.KPPg_TH.fillna(0), "Z")
    / metrics.cellvol
)

shfluxminusswrad = (
    -1
    * (surf.TFLUX - surf.oceQsw)
    / (rhoConst * Cp * metrics.dRF.isel(RC=0) * metrics.hFacC.isel(RC=0))
)

SW = (
    -1 * surf.oceQsw / (rhoConst * Cp) / (metrics.dRF[0] * metrics.hFacC[0]) * swprofile
)

# residual is LHS - RHS
# this is non-zero because of some global correction at the surface
budgets["RESIDUAL_TH"] = -1 * expand_depth(
    (
        shfluxminusswrad
        + (budgets.surf_corr_tend + budgets.MIX_TH + budgets.ADV_TH + SW).isel(RC=0)
    )
    - budgets.TOTTTEND.isel(RC=0) / 86400
)

budgets["SSHx"] = -9.81 * grid.diff(ssh.ETAN, "X") / metrics.DXC.data
budgets["SSHy"] = -9.81 * grid.diff(ssh.ETAN, "Y") / metrics.DYC.data

budgets["Um_Impl"] = (
    grid.diff(budgets.VISrI_Um.fillna(0), "Z")
    / metrics.rAw
    / metrics.drF
    / metrics.hFacW
)
budgets["Vm_Impl"] = (
    grid.diff(budgets.VISrI_Vm.fillna(0), "Z")
    / metrics.rAs
    / metrics.drF
    / metrics.hFacS
)

# Has LHS sign because we want to subtract it from the Advective part
budgets["Um_cor"] = -1 * (+1 * f0 * grid.interp(state.v, ("X", "Y")))
budgets["Vm_cor"] = -1 * (-1 * f0 * grid.interp(state.u, ("X", "Y")))

budgets = budgets.cf.guess_coord_axis()
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.8/site-packages/dask/array/core.py:377: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison
  o = func(*args, **kwargs)
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.8/site-packages/dask/array/core.py:377: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison
  o = func(*args, **kwargs)

Tendency terms#

Now subset the RHS terms to domain; add sum to total RHS tendencies.

Stupidly, the temperature flux divergences in budgets are with LHS sign; so invert that at the end

subrhs = subsetter(budgets)
rhs = xr.Dataset()
# these are names of variables in budgets that are added together in the for loop
var_terms = {
    "u": ["Um_dPHdx", "Um_Advec", "Um_Diss", "AB_gU", "SSHx", "Um_cor"],
    "v": ["Vm_dPHdy", "Vm_Advec", "Vm_Diss", "AB_gV", "SSHy", "Vm_cor"],
    "temp": ["ADV_TH", "surf_corr_tend", "RESIDUAL_TH"],
}

for var, varnames in var_terms.items():
    rhs[var] = (
        subrhs[varnames].pipe(interpolate, newz).to_array("term").sum("term")
    ) * 86400
    rhs[var].attrs["description"] = f"sum over {varnames}"
rhs
<xarray.Dataset>
Dimensions:  (RC: 216, time: 28)
Coordinates:
    XG       float64 -140.0
    XC       float64 -140.0
    YG       float64 3.607e-14
    YC       float64 0.025
  * time     (time) datetime64[ns] 2008-10-24T12:00:00 ... 2008-11-20T12:00:00
  * RC       (RC) float64 -107.8 -107.2 -106.8 -106.2 ... -1.25 -0.75 -0.25
Data variables:
    u        (time, RC) float64 -0.0001335 -0.0001465 ... 0.07101 0.07279
    v        (time, RC) float64 -0.0348 -0.03506 -0.03533 ... 0.007359 0.006686
    temp     (time, RC) float64 0.1456 0.1421 0.1386 ... -0.04826 -0.05021
write_to_txt(
    rhs.rename(renamer).compute().sel(time=slice(t0, None)),
    outdir=outdir,
    prefix="ocean_colpsh_",
    interleave=True,
    t0=t0,
)
_images/87cd382f55db3508c41a487d07ad804fa37544400c30f5ac2c747976e66182a7.png

initial condition#

subds = (
    state[["u", "v", "theta"]]
    .rename({"theta": "temp"})
    .pipe(subsetter)
    .pipe(interpolate, newz)
    .rename(renamer)
)
substate = (subds[["U", "W", "T"]].compute()) - offset

ic = substate.cf.interp(time=pd.Timestamp(t0))
write_to_txt(ic.compute(), outdir=outdir, prefix="ocean_ic_", interleave=False)
_images/9dda6a77157504960dc513d77169193596695a393a610813299b06edd508686d.png

restoring#

write_to_txt(substate, outdir=outdir, prefix="ocean_colfrc_", interleave=True, t0=t0)
_images/50a7781810589fb995311e42d6289715a89b366310be622fd60bbcdc86cfb534.png

Estimate surface forcing from JRA using coare3.6#

jrafull = pump.obs.read_jra_20()
# undo the local time correction for LES stuff
jrafull["time"] = jrafull.time + pd.Timedelta("7h")

# interpolate to location + subset in time
jrai = (
    jrafull.sel(
        time=slice(
            t0 - pd.Timedelta("2D"),  # something is buggy so start an extra day early
            substate.time[-1] + pd.Timedelta("23H"),
        )
    )
    .interp(latitude=substate.YC.data, longitude=substate.XC.data)
    .compute()
    .interpolate_na("time")
)
if jrai.sizes["time"] == 0:
    raise ValueError("JRA Dataset does not overlap with simulation time period!")

# calculate coare3.6 fluxes
coare, _ = pump.coare_fluxes_jra(
    subsetter(state, time=False).drop(["XG", "YG"]).interp(time=jrai.time.data),
    forcing=jrai,
)
coare.load()
daily_coare = coare.resample(time="D", loffset="12H").mean()

Estimate daily residuals

flux_residuals = xr.Dataset()
flux_residuals["short"] = daily_coare.short - surf["oceQsw"].pipe(subsetter, time=False)
flux_residuals["net"] = daily_coare.netflux - surf["TFLUX"].pipe(subsetter, time=False)
flux_residuals["taux"] = daily_coare.taux - surf.taux.pipe(subsetter, time=False)
flux_residuals["tauy"] = daily_coare.tauy - surf.tauy.pipe(subsetter, time=False)
flux_residuals = flux_residuals.reindex(
    time=coare.time.sel(
        time=slice(t0 - pd.Timedelta("1D"), t0 + pd.Timedelta(f"{les_time_length-3}D"))
    ),
    method="nearest",
).load()
for k, v in flux_residuals.isnull().sum().compute().items():
    assert v == 0, "NaNs in flux_residuals"

f = pump.les.debug_plots(flux_residuals)
f.savefig(f"{outdir}/coare_flux_residuals.png", dpi=150)
_images/eaba2860ba7ced8550b968cb08e5e7c9aff2bfba95e81f1ca07a7ef3a8bbdbcc.png
frc = xr.Dataset()
frc["sustr"] = coare.taux - flux_residuals.taux  # ρ * Δz * external forcing term
frc["svstr"] = coare.tauy - flux_residuals.tauy
frc["swrad"] = coare.short - flux_residuals.short
frc["shfluxminusswrad"] = (coare.netflux - coare.short) - (
    flux_residuals.net - flux_residuals.short
)

frc = frc.sel(time=slice(t0 - pd.Timedelta("1.5H"), None))
frc = frc.rename({"time": "time_index"})

time = frc.cf.indexes["time_index"].to_numpy()
frc["time"] = ("time_index", (time - t0).astype("timedelta64[s]").astype(float))
frc["time_index"].attrs["axis"] = "T"
frc.load()
write_to_txt(frc, outdir=outdir, prefix="ocean_frc_", interleave=False)
_images/f7afeeb666645fad8b4509de8c9586cb9cba73f505035389890106a422f16350.png

Test: coare fluxes vs mitgcm#

(surf["oceQsw"].pipe(subsetter)).plot()
daily_coare.short.plot()
(frc["swrad"].cf.resample(T="D", loffset="12H").mean()[:-1]).plot(marker="o", ls="none")
plt.legend(["MITgcm output", "COARE 3.6 daily", "corrected (ocean_frc*.txt)"])
plt.xlim((t0, None))
(14176.5, 14206.0)
_images/ed8d763101389d5e3b3b68c7429921aaaf8ed8546bdb182211d4457aa96f71f6.png
surf.taux.pipe(subsetter).plot()
daily_coare.taux.plot()
frc.cf.resample(T="D", loffset="12H").mean().sustr.plot(marker="o", ls="none")
plt.legend(["MITgcm output", "COARE 3.6 daily", "corrected (ocean_frc*.txt)"])
plt.xlim((t0, None))
(14176.5, 14206.0)
_images/95d76bb0ec19efbbb57ca93bb014953a9c1d3aa9d7c394998dcaffc21983d635.png
surf.tauy.pipe(subsetter).plot()
daily_coare.tauy.plot()
frc.cf.resample(T="D", loffset="12H").mean().svstr.plot(marker="o", ls="none")
plt.legend(["MITgcm output", "COARE 3.6 daily", "corrected (ocean_frc*.txt)"])
<matplotlib.legend.Legend at 0x2b56747a20a0>
_images/70e6cdaeead6dc6b4e9b6ec0fbed611f23c3c2f8e0ca3c27f26102e3fc1e6479.png
(surf.TFLUX - surf.oceQsw).pipe(subsetter).plot()
(daily_coare.netflux - daily_coare.short).plot()
frc.shfluxminusswrad.cf.resample(T="D", loffset="12H").mean().plot(
    marker="o", ls="none"
)
plt.legend(["MITgcm output", "COARE 3.6 daily", "corrected (ocean_frc*.txt)"])
<matplotlib.legend.Legend at 0x2b5674a304f0>
_images/a88d0f27834fd24e5da27575fd4957ba5732a9b0aebf3fd0e419e387718b148c.png

Raise error here so you don’t go past if running every cell :)

stop here
  File "<ipython-input-14-a96ba3aab008>", line 1
    stop here
         ^
SyntaxError: invalid syntax


Debugging checks#

Test: heat budget#

LHS = budgets.TOTTTEND / 86400
RHS = (
    expand_depth(shfluxminusswrad)
    + budgets.surf_corr_tend
    + budgets.MIX_TH
    + budgets.ADV_TH
    + SW
    + budgets.RESIDUAL_TH
)

# t =0 heat budget is weird at the surface
(LHS - RHS).isel(time=1, RC=[0, 10]).plot(robust=True, col="RC")
/glade/u/home/dcherian/miniconda3/envs/dcpy/lib/python3.8/site-packages/dask/array/core.py:377: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison
  o = func(*args, **kwargs)
<xarray.plot.facetgrid.FacetGrid at 0x2abcfc9f3610>
_images/1846199e6292517af9495f76b8fb24e1628497eed07d481f8f9a45049feb2a9e.png
residual = (LHS - RHS).isel(RC=0)
resmedian = residual.chunk({"time": 1}).median(["XC", "YC"]).compute()
rest = residual.mean(["XC", "YC"]).compute()
import matplotlib.pyplot as plt

plt.figure(constrained_layout=True)
rest.plot()
resmedian.plot()
plt.legend(["mean", "median"])
plt.title("heat budget residual at surface | from global correction")
plt.savefig("../images/heat-budget-residual.png")
_images/b1431bfc66786f82ae59cb971301a8b7263600ea3e72f8aa6d565c0db0459669.png
residual.isel(time=[1, 10, 20, 30]).plot(col="time", robust=True)
<xarray.plot.facetgrid.FacetGrid at 0x2b41fcb93970>
_images/9058fc027f156a7e12e82e25595afe2bae9a3fb1936edac12c0226cb53b9a987.png

Test : U budget#

RHS = (
    budgets.SSHx
    + budgets.Um_dPHdx
    + budgets.Um_Advec
    # + budget.Um_Cori  # is in Um_Advec
    + budgets.Um_Diss
    + budgets.Um_Impl
    + budgets.Um_Ext.fillna(0)
    + budgets.AB_gU
)

LHS = budgets.TOTUTEND / 86400
(LHS - RHS).isel(time=0, RC=[0, 10]).plot(robust=True, col="RC")
<xarray.plot.facetgrid.FacetGrid at 0x2b4291fe68b0>
_images/bcdae8581843210c4cd74420daebe966c869f96ba7601ea746ebc92ed1c3e835.png

test: V budget#

RHS = (
    budgets.SSHy
    + budgets.Vm_dPHdy
    + budgets.Vm_Advec
    # + budget.Um_Cori  # is in Um_Advec
    + budgets.Vm_Diss
    + budgets.Vm_Impl
    + budgets.Vm_Ext.fillna(0)
    + budgets.AB_gV
)

LHS = budgets.TOTVTEND / 86400
(LHS - RHS).isel(time=0, RC=[0, 10]).plot(robust=True, col="RC")
<xarray.plot.facetgrid.FacetGrid at 0x2b429365b220>
_images/0314ceaed9f1aa45f9c4ede9cb232af797f21f466195ca6a6dc6b8117b91de21.png

ROMS#

Todo /glade/p/cgd/oce/people/dwhitt/TPOS/tpos_LES_runs_setup_scripts/tpos-DIABLO/diablo_2.0/tpos_config_files/ROMS_PSH_6HRLIN_0N140W_108x108x288_5OCT2020_fixedeps/

  1. Staggered grid

  2. Update grid

  3. something is wrong with your ICs (ocean_ic_S.txt and ocean_ic_T.txt and ocean_ic_U.txt). The values are much too small, except for U = -1 + epsilon. U should be m/s (anomaly). S psu (anomaly). T deg C or Kelvin (anomaly).

  4. except ocean_frc_time.txt. Something is off there with the units/order of magnitude. You need to make sure the time that you want the LES to start is zero, so you have to subtract off the initial time and then have all times as seconds after (and before) that point. I think it is necessary to have one time index (i.e. one negative ocean_frc_time value) before the targeted LES start time for interpolation.

  5. please add the time step before each profile as in ocean_colpsh_T.txt in this directory:/glade/p/cgd/oce/people/dwhitt/TPOS/tpos_LES_runs_setup_scripts/tpos-DIABLO/diablo_2.0/tpos_config_files

  6. put the outputs in units of deg C per day, m/s per day, psu per day. (I changed the LES to take in units of per day to mitigate precision issues).

  7. You need ocean_colfrc_*.txt files, which are for restoring. They are just the basic fields, T(z,t), U(z,t), S(z,t) and V(z,t) (called “W” in LES world; don’t ask…).

import xroms

outdir = "/glade/work/dcherian/pump/les_forcing/roms_test/"
dirname = "/glade/p/cgd/oce/people/dwhitt/TPOS/tpos20/OUT/"
files = [
    f"{dirname}/ocean_dia_tpos20_3_d2_2.nc",
    f"{dirname}/ocean_avg_tpos20_3_d2_2.nc",
]

chunk = {"xi": 500, "eta": 120}

ds0 = pump.model.model.read_roms_dataset(files, xi=20, eta=500, ocean_time=50)
rds, grid = xroms.roms_dataset(ds0)
rds["u"] = grid.interp(rds.u, axis="X", boundary="extend")
rds["v"] = grid.interp(rds.v, axis="Y", boundary="extend")
romsds = pump.model.model.make_cartesian(rds)
romsds
/glade/u/home/dcherian/pump/pump/model/model.py:224: RuntimeWarning: Converting a CFTimeIndex with dates from a non-standard calendar, 'noleap', to a pandas.DatetimeIndex, which uses dates from the standard calendar.  This may lead to subtle errors in operations that depend on the length of time between dates.
  ds["ocean_time"] = ds.indexes["ocean_time"].to_datetimeindex()
<xarray.Dataset>
Dimensions:        (boundary: 4, lat_rho: 481, lat_v: 480, lon_rho: 1501, lon_u: 1500, ocean_time: 140, s_rho: 50, s_w: 51, tracer: 2)
Coordinates:
    ntimes         int32 25201
    ndtfast        int32 20
    dt             float64 120.0
    dtfast         float64 6.0
    dstart         object 0027-09-02 00:00:00
    nHIS           int32 3600
    ndefHIS        int32 0
    nRST           int32 3600
    ntsAVG         int32 1
    nAVG           int32 180
    ndefAVG        int32 0
    ntsDIA         int32 1
    nDIA           int32 180
    ndefDIA        int32 0
    nSTA           int32 30
    Falpha         float64 2.0
    Fgamma         float64 0.284
    nl_tnu4        (tracer) float64 dask.array<chunksize=(2,), meta=np.ndarray>
    nl_visc4       float64 1.24e+09
    LuvSponge      int32 0
    LtracerSponge  (tracer) int32 dask.array<chunksize=(2,), meta=np.ndarray>
    Akt_bak        (tracer) float64 dask.array<chunksize=(2,), meta=np.ndarray>
    Akv_bak        float64 1e-05
    rdrg           float64 0.000265
    rdrg2          float64 0.003
    Zob            float64 0.02
    Zos            float64 0.02
    Znudg          float64 0.0
    M2nudg         float64 0.0
    M3nudg         float64 1.0
    Tnudg          (tracer) float64 dask.array<chunksize=(2,), meta=np.ndarray>
    FSobc_in       (boundary) float64 dask.array<chunksize=(4,), meta=np.ndarray>
    FSobc_out      (boundary) float64 dask.array<chunksize=(4,), meta=np.ndarray>
    M2obc_in       (boundary) float64 dask.array<chunksize=(4,), meta=np.ndarray>
    M2obc_out      (boundary) float64 dask.array<chunksize=(4,), meta=np.ndarray>
    Tobc_in        (boundary, tracer) float64 dask.array<chunksize=(4, 2), meta=np.ndarray>
    Tobc_out       (boundary, tracer) float64 dask.array<chunksize=(4, 2), meta=np.ndarray>
    M3obc_in       (boundary) float64 dask.array<chunksize=(4,), meta=np.ndarray>
    M3obc_out      (boundary) float64 dask.array<chunksize=(4,), meta=np.ndarray>
    rho0           float64 1.025e+03
    gamma2         float64 1.0
    LuvSrc         int32 0
    LwSrc          int32 0
    LtracerSrc     (tracer) int32 dask.array<chunksize=(2,), meta=np.ndarray>
    LsshCLM        int32 0
    Lm2CLM         int32 0
    Lm3CLM         int32 0
    LtracerCLM     (tracer) int32 dask.array<chunksize=(2,), meta=np.ndarray>
    LnudgeM2CLM    int32 0
    LnudgeM3CLM    int32 0
    LnudgeTCLM     (tracer) int32 dask.array<chunksize=(2,), meta=np.ndarray>
    Lm2PSH         int32 0
    Lm3PSH         int32 0
    LtracerPSH     (tracer) int32 dask.array<chunksize=(2,), meta=np.ndarray>
    spherical      int32 1
    xl             float64 8.163e+06
    el             float64 2.674e+06
    Vtransform     int32 1
    Vstretching    int32 1
    Tcline         float64 75.0
    hc             float64 75.0
    grid           int32 1
  * s_rho          (s_rho) float64 -0.99 -0.97 -0.95 -0.93 ... -0.05 -0.03 -0.01
  * s_w            (s_w) float64 -1.0 -0.98 -0.96 -0.94 ... -0.04 -0.02 0.0
    Cs_r           (s_rho) float64 dask.array<chunksize=(50,), meta=np.ndarray>
    Cs_w           (s_w) float64 dask.array<chunksize=(51,), meta=np.ndarray>
    h              (lat_rho, lon_rho) float64 dask.array<chunksize=(481, 20), meta=np.ndarray>
    f              (lat_rho, lon_rho) float64 dask.array<chunksize=(481, 20), meta=np.ndarray>
    pm             (lat_rho, lon_rho) float64 dask.array<chunksize=(481, 20), meta=np.ndarray>
    pn             (lat_rho, lon_rho) float64 dask.array<chunksize=(481, 20), meta=np.ndarray>
  * lon_rho        (lon_rho) float64 -170.0 -170.0 -169.9 ... -95.1 -95.05 -95.0
  * lat_rho        (lat_rho) float64 -12.0 -11.95 -11.9 ... 11.9 11.95 12.0
  * lon_u          (lon_u) float64 -170.0 -169.9 -169.9 ... -95.12 -95.07 -95.02
    lat_u          (lat_rho, lon_u) float64 dask.array<chunksize=(481, 20), meta=np.ndarray>
    lon_v          (lat_v, lon_rho) float64 dask.array<chunksize=(480, 20), meta=np.ndarray>
  * lat_v          (lat_v) float64 -11.97 -11.92 -11.87 ... 11.87 11.92 11.97
    lon_psi        (lat_v, lon_u) float64 dask.array<chunksize=(480, 20), meta=np.ndarray>
    lat_psi        (lat_v, lon_u) float64 dask.array<chunksize=(480, 20), meta=np.ndarray>
    angle          (lat_rho, lon_rho) float64 dask.array<chunksize=(481, 20), meta=np.ndarray>
    mask_rho       (lat_rho, lon_rho) float64 dask.array<chunksize=(481, 20), meta=np.ndarray>
    mask_u         (lat_rho, lon_u) float64 dask.array<chunksize=(481, 20), meta=np.ndarray>
    mask_v         (lat_v, lon_rho) float64 dask.array<chunksize=(480, 20), meta=np.ndarray>
    mask_psi       (lat_v, lon_u) float64 dask.array<chunksize=(480, 20), meta=np.ndarray>
  * ocean_time     (ocean_time) datetime64[ns] 1985-10-02T03:00:00 ... 1985-1...
    z_w            (ocean_time, s_w, lat_rho, lon_rho) float64 dask.array<chunksize=(50, 51, 481, 20), meta=np.ndarray>
    z_w_u          (ocean_time, s_w, lat_rho, lon_u) float64 dask.array<chunksize=(50, 51, 481, 19), meta=np.ndarray>
    z_w_v          (ocean_time, s_w, lat_v, lon_rho) float64 dask.array<chunksize=(50, 51, 480, 20), meta=np.ndarray>
    z_w_psi        (ocean_time, s_w, lat_v, lon_u) float64 dask.array<chunksize=(50, 51, 480, 19), meta=np.ndarray>
    z_rho          (ocean_time, s_rho, lat_rho, lon_rho) float64 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    z_rho_u        (ocean_time, s_rho, lat_rho, lon_u) float64 dask.array<chunksize=(50, 50, 481, 19), meta=np.ndarray>
    z_rho_v        (ocean_time, s_rho, lat_v, lon_rho) float64 dask.array<chunksize=(50, 50, 480, 20), meta=np.ndarray>
    z_rho_psi      (ocean_time, s_rho, lat_v, lon_u) float64 dask.array<chunksize=(50, 50, 480, 19), meta=np.ndarray>
    z_rho0         (s_rho, lat_rho, lon_rho) float64 dask.array<chunksize=(50, 481, 20), meta=np.ndarray>
    z_rho_u0       (s_rho, lat_rho, lon_u) float64 dask.array<chunksize=(50, 481, 19), meta=np.ndarray>
    z_rho_v0       (s_rho, lat_v, lon_rho) float64 dask.array<chunksize=(50, 480, 20), meta=np.ndarray>
    z_rho_psi0     (s_rho, lat_v, lon_u) float64 dask.array<chunksize=(50, 480, 19), meta=np.ndarray>
    z_w0           (s_w, lat_rho, lon_rho) float64 dask.array<chunksize=(51, 481, 20), meta=np.ndarray>
    z_w_u0         (s_w, lat_rho, lon_u) float64 dask.array<chunksize=(51, 481, 19), meta=np.ndarray>
    z_w_v0         (s_w, lat_v, lon_rho) float64 dask.array<chunksize=(51, 480, 20), meta=np.ndarray>
    z_w_psi0       (s_w, lat_v, lon_u) float64 dask.array<chunksize=(51, 480, 19), meta=np.ndarray>
Dimensions without coordinates: boundary, tracer
Data variables:
    zeta           (ocean_time, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 481, 20), meta=np.ndarray>
    ubar           (ocean_time, lat_rho, lon_u) float32 dask.array<chunksize=(50, 481, 20), meta=np.ndarray>
    vbar           (ocean_time, lat_v, lon_rho) float32 dask.array<chunksize=(50, 480, 20), meta=np.ndarray>
    u              (ocean_time, s_rho, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 50, 481, 1), meta=np.ndarray>
    v              (ocean_time, s_rho, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 50, 1, 20), meta=np.ndarray>
    w              (ocean_time, s_w, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 51, 481, 20), meta=np.ndarray>
    temp           (ocean_time, s_rho, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    salt           (ocean_time, s_rho, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    Hsbl           (ocean_time, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 481, 20), meta=np.ndarray>
    AKv            (ocean_time, s_w, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 51, 481, 20), meta=np.ndarray>
    AKt            (ocean_time, s_w, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 51, 481, 20), meta=np.ndarray>
    shflux         (ocean_time, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 481, 20), meta=np.ndarray>
    ssflux         (ocean_time, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 481, 20), meta=np.ndarray>
    latent         (ocean_time, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 481, 20), meta=np.ndarray>
    sensible       (ocean_time, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 481, 20), meta=np.ndarray>
    lwrad          (ocean_time, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 481, 20), meta=np.ndarray>
    evaporation    (ocean_time, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 481, 20), meta=np.ndarray>
    rain           (ocean_time, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 481, 20), meta=np.ndarray>
    swrad          (ocean_time, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 481, 20), meta=np.ndarray>
    sustr          (ocean_time, lat_rho, lon_u) float32 dask.array<chunksize=(50, 481, 20), meta=np.ndarray>
    svstr          (ocean_time, lat_v, lon_rho) float32 dask.array<chunksize=(50, 480, 20), meta=np.ndarray>
    u_cor          (ocean_time, s_rho, lat_rho, lon_u) float32 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    v_cor          (ocean_time, s_rho, lat_v, lon_rho) float32 dask.array<chunksize=(50, 50, 480, 20), meta=np.ndarray>
    u_vadv         (ocean_time, s_rho, lat_rho, lon_u) float32 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    v_vadv         (ocean_time, s_rho, lat_v, lon_rho) float32 dask.array<chunksize=(50, 50, 480, 20), meta=np.ndarray>
    u_hadv         (ocean_time, s_rho, lat_rho, lon_u) float32 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    v_hadv         (ocean_time, s_rho, lat_v, lon_rho) float32 dask.array<chunksize=(50, 50, 480, 20), meta=np.ndarray>
    u_xadv         (ocean_time, s_rho, lat_rho, lon_u) float32 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    v_xadv         (ocean_time, s_rho, lat_v, lon_rho) float32 dask.array<chunksize=(50, 50, 480, 20), meta=np.ndarray>
    u_yadv         (ocean_time, s_rho, lat_rho, lon_u) float32 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    v_yadv         (ocean_time, s_rho, lat_v, lon_rho) float32 dask.array<chunksize=(50, 50, 480, 20), meta=np.ndarray>
    u_prsgrd       (ocean_time, s_rho, lat_rho, lon_u) float32 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    v_prsgrd       (ocean_time, s_rho, lat_v, lon_rho) float32 dask.array<chunksize=(50, 50, 480, 20), meta=np.ndarray>
    u_vvisc        (ocean_time, s_rho, lat_rho, lon_u) float32 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    v_vvisc        (ocean_time, s_rho, lat_v, lon_rho) float32 dask.array<chunksize=(50, 50, 480, 20), meta=np.ndarray>
    u_hvisc        (ocean_time, s_rho, lat_rho, lon_u) float32 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    v_hvisc        (ocean_time, s_rho, lat_v, lon_rho) float32 dask.array<chunksize=(50, 50, 480, 20), meta=np.ndarray>
    u_xvisc        (ocean_time, s_rho, lat_rho, lon_u) float32 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    v_xvisc        (ocean_time, s_rho, lat_v, lon_rho) float32 dask.array<chunksize=(50, 50, 480, 20), meta=np.ndarray>
    u_yvisc        (ocean_time, s_rho, lat_rho, lon_u) float32 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    v_yvisc        (ocean_time, s_rho, lat_v, lon_rho) float32 dask.array<chunksize=(50, 50, 480, 20), meta=np.ndarray>
    u_accel        (ocean_time, s_rho, lat_rho, lon_u) float32 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    v_accel        (ocean_time, s_rho, lat_v, lon_rho) float32 dask.array<chunksize=(50, 50, 480, 20), meta=np.ndarray>
    temp_hadv      (ocean_time, s_rho, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    temp_xadv      (ocean_time, s_rho, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    temp_yadv      (ocean_time, s_rho, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    temp_vadv      (ocean_time, s_rho, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    temp_hdiff     (ocean_time, s_rho, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    temp_xdiff     (ocean_time, s_rho, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    temp_ydiff     (ocean_time, s_rho, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    temp_sdiff     (ocean_time, s_rho, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    temp_vdiff     (ocean_time, s_rho, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    temp_rate      (ocean_time, s_rho, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    salt_hadv      (ocean_time, s_rho, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    salt_xadv      (ocean_time, s_rho, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    salt_yadv      (ocean_time, s_rho, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    salt_vadv      (ocean_time, s_rho, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    salt_hdiff     (ocean_time, s_rho, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    salt_xdiff     (ocean_time, s_rho, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    salt_ydiff     (ocean_time, s_rho, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    salt_sdiff     (ocean_time, s_rho, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    salt_vdiff     (ocean_time, s_rho, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    salt_rate      (ocean_time, s_rho, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    pm_v           (lat_v, lon_rho) float64 dask.array<chunksize=(480, 20), meta=np.ndarray>
    pn_u           (lat_rho, lon_u) float64 dask.array<chunksize=(481, 19), meta=np.ndarray>
    pm_u           (lat_rho, lon_u) float64 dask.array<chunksize=(481, 19), meta=np.ndarray>
    pn_v           (lat_v, lon_rho) float64 dask.array<chunksize=(480, 20), meta=np.ndarray>
    pm_psi         (lat_v, lon_u) float64 dask.array<chunksize=(480, 19), meta=np.ndarray>
    pn_psi         (lat_v, lon_u) float64 dask.array<chunksize=(480, 19), meta=np.ndarray>
    dx             (lat_rho, lon_rho) float64 dask.array<chunksize=(481, 20), meta=np.ndarray>
    dx_u           (lat_rho, lon_u) float64 dask.array<chunksize=(481, 19), meta=np.ndarray>
    dx_v           (lat_v, lon_rho) float64 dask.array<chunksize=(480, 20), meta=np.ndarray>
    dx_psi         (lat_v, lon_u) float64 dask.array<chunksize=(480, 19), meta=np.ndarray>
    dy             (lat_rho, lon_rho) float64 dask.array<chunksize=(481, 20), meta=np.ndarray>
    dy_u           (lat_rho, lon_u) float64 dask.array<chunksize=(481, 19), meta=np.ndarray>
    dy_v           (lat_v, lon_rho) float64 dask.array<chunksize=(480, 20), meta=np.ndarray>
    dy_psi         (lat_v, lon_u) float64 dask.array<chunksize=(480, 19), meta=np.ndarray>
    dz             (ocean_time, s_rho, lat_rho, lon_rho) float64 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    dz_w           (ocean_time, s_w, lat_rho, lon_rho) float64 dask.array<chunksize=(50, 1, 481, 20), meta=np.ndarray>
    dz_u           (ocean_time, s_rho, lat_rho, lon_u) float64 dask.array<chunksize=(50, 50, 481, 19), meta=np.ndarray>
    dz_w_u         (ocean_time, s_w, lat_rho, lon_u) float64 dask.array<chunksize=(50, 1, 481, 19), meta=np.ndarray>
    dz_v           (ocean_time, s_rho, lat_v, lon_rho) float64 dask.array<chunksize=(50, 50, 480, 20), meta=np.ndarray>
    dz_w_v         (ocean_time, s_w, lat_v, lon_rho) float64 dask.array<chunksize=(50, 1, 480, 20), meta=np.ndarray>
    dz_psi         (ocean_time, s_rho, lat_v, lon_u) float64 dask.array<chunksize=(50, 50, 480, 19), meta=np.ndarray>
    dz_w_psi       (ocean_time, s_w, lat_v, lon_u) float64 dask.array<chunksize=(50, 1, 480, 19), meta=np.ndarray>
    dz0            (s_rho, lat_rho, lon_rho) float64 dask.array<chunksize=(50, 481, 20), meta=np.ndarray>
    dz_w0          (s_w, lat_rho, lon_rho) float64 dask.array<chunksize=(1, 481, 20), meta=np.ndarray>
    dz_u0          (s_rho, lat_rho, lon_u) float64 dask.array<chunksize=(50, 481, 19), meta=np.ndarray>
    dz_w_u0        (s_w, lat_rho, lon_u) float64 dask.array<chunksize=(1, 481, 19), meta=np.ndarray>
    dz_v0          (s_rho, lat_v, lon_rho) float64 dask.array<chunksize=(50, 480, 20), meta=np.ndarray>
    dz_w_v0        (s_w, lat_v, lon_rho) float64 dask.array<chunksize=(1, 480, 20), meta=np.ndarray>
    dz_psi0        (s_rho, lat_v, lon_u) float64 dask.array<chunksize=(50, 480, 19), meta=np.ndarray>
    dz_w_psi0      (s_w, lat_v, lon_u) float64 dask.array<chunksize=(1, 480, 19), meta=np.ndarray>
    dA             (lat_rho, lon_rho) float64 dask.array<chunksize=(481, 20), meta=np.ndarray>
    dA_u           (lat_rho, lon_u) float64 dask.array<chunksize=(481, 19), meta=np.ndarray>
    dA_v           (lat_v, lon_rho) float64 dask.array<chunksize=(480, 20), meta=np.ndarray>
    dA_psi         (lat_v, lon_u) float64 dask.array<chunksize=(480, 19), meta=np.ndarray>
    dV             (ocean_time, s_rho, lat_rho, lon_rho) float64 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    dV_w           (ocean_time, s_w, lat_rho, lon_rho) float64 dask.array<chunksize=(50, 1, 481, 20), meta=np.ndarray>
    dV_u           (ocean_time, s_rho, lat_rho, lon_u) float64 dask.array<chunksize=(50, 50, 481, 19), meta=np.ndarray>
    dV_w_u         (ocean_time, s_w, lat_rho, lon_u) float64 dask.array<chunksize=(50, 1, 481, 19), meta=np.ndarray>
    dV_v           (ocean_time, s_rho, lat_v, lon_rho) float64 dask.array<chunksize=(50, 50, 480, 20), meta=np.ndarray>
    dV_w_v         (ocean_time, s_w, lat_v, lon_rho) float64 dask.array<chunksize=(50, 1, 480, 20), meta=np.ndarray>
    dV_psi         (ocean_time, s_rho, lat_v, lon_u) float64 dask.array<chunksize=(50, 50, 480, 19), meta=np.ndarray>
    dV_w_psi       (ocean_time, s_w, lat_v, lon_u) float64 dask.array<chunksize=(50, 1, 480, 19), meta=np.ndarray>
Attributes:
    file:              /glade/scratch/dwhitt/tpos20/OUT/ocean_dia_tpos20_3_d2...
    format:            netCDF-3 64bit offset file
    Conventions:       CF-1.4, SGRID-0.3
    type:              ROMS/TOMS diagnostics file
    title:             TPOS 1/20deg
    var_info:          /glade/p/nsc/ncgd0043/tpos_ROMS_runs_setup_scripts/tpo...
    rst_file:          /glade/scratch/dwhitt/tpos20/OUT/ocean_rst_tpos20_3_d2...
    avg_file:          /glade/scratch/dwhitt/tpos20/OUT/ocean_avg_tpos20_3_d2...
    dia_file:          /glade/scratch/dwhitt/tpos20/OUT/ocean_dia_tpos20_3_d2...
    sta_file:          /glade/scratch/dwhitt/tpos20/OUT/ocean_sta_tpos20_3_d2...
    grd_file:          /glade/p/nsc/ncgd0043/tpos_ROMS_runs_setup_scripts/tpo...
    ini_file:          /glade/p/nsc/ncgd0043/tpos20/OUT/ocean_rst_tpos20_3_d2.nc
    frc_file_01:       /glade/p/nsc/ncgd0043/tpos_ROMS_runs_setup_scripts/tpo...
    frc_file_02:       /glade/p/nsc/ncgd0043/tpos_ROMS_runs_setup_scripts/tpo...
    frc_file_03:       /glade/p/nsc/ncgd0043/tpos_ROMS_runs_setup_scripts/tpo...
    bry_file:          /glade/p/nsc/ncgd0043/tpos_ROMS_runs_setup_scripts/tpo...
    script_file:       ./tpos20_3_diag_2.in
    spos_file:         /glade/p/nsc/ncgd0043/tpos_ROMS_runs_setup_scripts/tpo...
    NLM_LBC:           \nEDGE:  WEST   SOUTH  EAST   NORTH  \nzeta:  Che    C...
    svn_url:           https:://myroms.org/svn/src
    svn_rev:           Unversioned directory
    code_dir:          /glade/p/nsc/ncgd0043/tpos_ROMS_runs_setup_scripts/ROM...
    header_dir:        /glade/p/nsc/ncgd0043/tpos_ROMS_runs_setup_scripts/tpo...
    header_file:       tpos20.h
    os:                Linux
    cpu:               x86_64
    compiler_system:   ifort
    compiler_command:  /glade/u/apps/ch/opt/ncarcompilers/0.5.0/intel/18.0.5/...
    compiler_flags:    -fp-model precise -ip -O3
    tiling:            018x012
    history:           ROMS/TOMS, Version 3.7, Tuesday - July 30, 2019 -  5:5...
    ana_file:          ROMS/Functionals/ana_btflux.h
    CPP_options:       TPOS20, ALBEDO_JRA, ANA_BSFLUX, ANA_BTFLUX, ASSUMED_SH...
romsds
<xarray.Dataset>
Dimensions:        (boundary: 4, lat_rho: 481, lat_v: 480, lon_rho: 1501, lon_u: 1500, ocean_time: 140, s_rho: 50, s_w: 51, tracer: 2)
Coordinates:
    ntimes         int32 25201
    ndtfast        int32 20
    dt             float64 120.0
    dtfast         float64 6.0
    dstart         object 0027-09-02 00:00:00
    nHIS           int32 3600
    ndefHIS        int32 0
    nRST           int32 3600
    ntsAVG         int32 1
    nAVG           int32 180
    ndefAVG        int32 0
    ntsDIA         int32 1
    nDIA           int32 180
    ndefDIA        int32 0
    nSTA           int32 30
    Falpha         float64 2.0
    Fgamma         float64 0.284
    nl_tnu4        (tracer) float64 dask.array<chunksize=(2,), meta=np.ndarray>
    nl_visc4       float64 1.24e+09
    LuvSponge      int32 0
    LtracerSponge  (tracer) int32 dask.array<chunksize=(2,), meta=np.ndarray>
    Akt_bak        (tracer) float64 dask.array<chunksize=(2,), meta=np.ndarray>
    Akv_bak        float64 1e-05
    rdrg           float64 0.000265
    rdrg2          float64 0.003
    Zob            float64 0.02
    Zos            float64 0.02
    Znudg          float64 0.0
    M2nudg         float64 0.0
    M3nudg         float64 1.0
    Tnudg          (tracer) float64 dask.array<chunksize=(2,), meta=np.ndarray>
    FSobc_in       (boundary) float64 dask.array<chunksize=(4,), meta=np.ndarray>
    FSobc_out      (boundary) float64 dask.array<chunksize=(4,), meta=np.ndarray>
    M2obc_in       (boundary) float64 dask.array<chunksize=(4,), meta=np.ndarray>
    M2obc_out      (boundary) float64 dask.array<chunksize=(4,), meta=np.ndarray>
    Tobc_in        (boundary, tracer) float64 dask.array<chunksize=(4, 2), meta=np.ndarray>
    Tobc_out       (boundary, tracer) float64 dask.array<chunksize=(4, 2), meta=np.ndarray>
    M3obc_in       (boundary) float64 dask.array<chunksize=(4,), meta=np.ndarray>
    M3obc_out      (boundary) float64 dask.array<chunksize=(4,), meta=np.ndarray>
    rho0           float64 1.025e+03
    gamma2         float64 1.0
    LuvSrc         int32 0
    LwSrc          int32 0
    LtracerSrc     (tracer) int32 dask.array<chunksize=(2,), meta=np.ndarray>
    LsshCLM        int32 0
    Lm2CLM         int32 0
    Lm3CLM         int32 0
    LtracerCLM     (tracer) int32 dask.array<chunksize=(2,), meta=np.ndarray>
    LnudgeM2CLM    int32 0
    LnudgeM3CLM    int32 0
    LnudgeTCLM     (tracer) int32 dask.array<chunksize=(2,), meta=np.ndarray>
    Lm2PSH         int32 0
    Lm3PSH         int32 0
    LtracerPSH     (tracer) int32 dask.array<chunksize=(2,), meta=np.ndarray>
    spherical      int32 1
    xl             float64 8.163e+06
    el             float64 2.674e+06
    Vtransform     int32 1
    Vstretching    int32 1
    Tcline         float64 75.0
    hc             float64 75.0
    grid           int32 1
  * s_rho          (s_rho) float64 -0.99 -0.97 -0.95 -0.93 ... -0.05 -0.03 -0.01
  * s_w            (s_w) float64 -1.0 -0.98 -0.96 -0.94 ... -0.04 -0.02 0.0
    Cs_r           (s_rho) float64 dask.array<chunksize=(50,), meta=np.ndarray>
    Cs_w           (s_w) float64 dask.array<chunksize=(51,), meta=np.ndarray>
    h              (lat_rho, lon_rho) float64 dask.array<chunksize=(481, 20), meta=np.ndarray>
    f              (lat_rho, lon_rho) float64 dask.array<chunksize=(481, 20), meta=np.ndarray>
    pm             (lat_rho, lon_rho) float64 dask.array<chunksize=(481, 20), meta=np.ndarray>
    pn             (lat_rho, lon_rho) float64 dask.array<chunksize=(481, 20), meta=np.ndarray>
  * lon_rho        (lon_rho) float64 -170.0 -170.0 -169.9 ... -95.1 -95.05 -95.0
  * lat_rho        (lat_rho) float64 -12.0 -11.95 -11.9 ... 11.9 11.95 12.0
  * lon_u          (lon_u) float64 -170.0 -169.9 -169.9 ... -95.12 -95.07 -95.02
    lat_u          (lat_rho, lon_u) float64 dask.array<chunksize=(481, 20), meta=np.ndarray>
    lon_v          (lat_v, lon_rho) float64 dask.array<chunksize=(480, 20), meta=np.ndarray>
  * lat_v          (lat_v) float64 -11.97 -11.92 -11.87 ... 11.87 11.92 11.97
    lon_psi        (lat_v, lon_u) float64 dask.array<chunksize=(480, 20), meta=np.ndarray>
    lat_psi        (lat_v, lon_u) float64 dask.array<chunksize=(480, 20), meta=np.ndarray>
    angle          (lat_rho, lon_rho) float64 dask.array<chunksize=(481, 20), meta=np.ndarray>
    mask_rho       (lat_rho, lon_rho) float64 dask.array<chunksize=(481, 20), meta=np.ndarray>
    mask_u         (lat_rho, lon_u) float64 dask.array<chunksize=(481, 20), meta=np.ndarray>
    mask_v         (lat_v, lon_rho) float64 dask.array<chunksize=(480, 20), meta=np.ndarray>
    mask_psi       (lat_v, lon_u) float64 dask.array<chunksize=(480, 20), meta=np.ndarray>
  * ocean_time     (ocean_time) datetime64[ns] 1985-10-02T03:00:00 ... 1985-1...
    z_w            (ocean_time, s_w, lat_rho, lon_rho) float64 dask.array<chunksize=(50, 51, 481, 20), meta=np.ndarray>
    z_w_u          (ocean_time, s_w, lat_rho, lon_u) float64 dask.array<chunksize=(50, 51, 481, 19), meta=np.ndarray>
    z_w_v          (ocean_time, s_w, lat_v, lon_rho) float64 dask.array<chunksize=(50, 51, 480, 20), meta=np.ndarray>
    z_w_psi        (ocean_time, s_w, lat_v, lon_u) float64 dask.array<chunksize=(50, 51, 480, 19), meta=np.ndarray>
    z_rho          (ocean_time, s_rho, lat_rho, lon_rho) float64 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    z_rho_u        (ocean_time, s_rho, lat_rho, lon_u) float64 dask.array<chunksize=(50, 50, 481, 19), meta=np.ndarray>
    z_rho_v        (ocean_time, s_rho, lat_v, lon_rho) float64 dask.array<chunksize=(50, 50, 480, 20), meta=np.ndarray>
    z_rho_psi      (ocean_time, s_rho, lat_v, lon_u) float64 dask.array<chunksize=(50, 50, 480, 19), meta=np.ndarray>
    z_rho0         (s_rho, lat_rho, lon_rho) float64 dask.array<chunksize=(50, 481, 20), meta=np.ndarray>
    z_rho_u0       (s_rho, lat_rho, lon_u) float64 dask.array<chunksize=(50, 481, 19), meta=np.ndarray>
    z_rho_v0       (s_rho, lat_v, lon_rho) float64 dask.array<chunksize=(50, 480, 20), meta=np.ndarray>
    z_rho_psi0     (s_rho, lat_v, lon_u) float64 dask.array<chunksize=(50, 480, 19), meta=np.ndarray>
    z_w0           (s_w, lat_rho, lon_rho) float64 dask.array<chunksize=(51, 481, 20), meta=np.ndarray>
    z_w_u0         (s_w, lat_rho, lon_u) float64 dask.array<chunksize=(51, 481, 19), meta=np.ndarray>
    z_w_v0         (s_w, lat_v, lon_rho) float64 dask.array<chunksize=(51, 480, 20), meta=np.ndarray>
    z_w_psi0       (s_w, lat_v, lon_u) float64 dask.array<chunksize=(51, 480, 19), meta=np.ndarray>
Dimensions without coordinates: boundary, tracer
Data variables:
    zeta           (ocean_time, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 481, 20), meta=np.ndarray>
    ubar           (ocean_time, lat_rho, lon_u) float32 dask.array<chunksize=(50, 481, 20), meta=np.ndarray>
    vbar           (ocean_time, lat_v, lon_rho) float32 dask.array<chunksize=(50, 480, 20), meta=np.ndarray>
    u              (ocean_time, s_rho, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 50, 481, 1), meta=np.ndarray>
    v              (ocean_time, s_rho, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 50, 1, 20), meta=np.ndarray>
    w              (ocean_time, s_w, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 51, 481, 20), meta=np.ndarray>
    temp           (ocean_time, s_rho, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    salt           (ocean_time, s_rho, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    Hsbl           (ocean_time, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 481, 20), meta=np.ndarray>
    AKv            (ocean_time, s_w, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 51, 481, 20), meta=np.ndarray>
    AKt            (ocean_time, s_w, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 51, 481, 20), meta=np.ndarray>
    shflux         (ocean_time, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 481, 20), meta=np.ndarray>
    ssflux         (ocean_time, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 481, 20), meta=np.ndarray>
    latent         (ocean_time, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 481, 20), meta=np.ndarray>
    sensible       (ocean_time, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 481, 20), meta=np.ndarray>
    lwrad          (ocean_time, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 481, 20), meta=np.ndarray>
    evaporation    (ocean_time, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 481, 20), meta=np.ndarray>
    rain           (ocean_time, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 481, 20), meta=np.ndarray>
    swrad          (ocean_time, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 481, 20), meta=np.ndarray>
    sustr          (ocean_time, lat_rho, lon_u) float32 dask.array<chunksize=(50, 481, 20), meta=np.ndarray>
    svstr          (ocean_time, lat_v, lon_rho) float32 dask.array<chunksize=(50, 480, 20), meta=np.ndarray>
    u_cor          (ocean_time, s_rho, lat_rho, lon_u) float32 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    v_cor          (ocean_time, s_rho, lat_v, lon_rho) float32 dask.array<chunksize=(50, 50, 480, 20), meta=np.ndarray>
    u_vadv         (ocean_time, s_rho, lat_rho, lon_u) float32 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    v_vadv         (ocean_time, s_rho, lat_v, lon_rho) float32 dask.array<chunksize=(50, 50, 480, 20), meta=np.ndarray>
    u_hadv         (ocean_time, s_rho, lat_rho, lon_u) float32 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    v_hadv         (ocean_time, s_rho, lat_v, lon_rho) float32 dask.array<chunksize=(50, 50, 480, 20), meta=np.ndarray>
    u_xadv         (ocean_time, s_rho, lat_rho, lon_u) float32 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    v_xadv         (ocean_time, s_rho, lat_v, lon_rho) float32 dask.array<chunksize=(50, 50, 480, 20), meta=np.ndarray>
    u_yadv         (ocean_time, s_rho, lat_rho, lon_u) float32 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    v_yadv         (ocean_time, s_rho, lat_v, lon_rho) float32 dask.array<chunksize=(50, 50, 480, 20), meta=np.ndarray>
    u_prsgrd       (ocean_time, s_rho, lat_rho, lon_u) float32 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    v_prsgrd       (ocean_time, s_rho, lat_v, lon_rho) float32 dask.array<chunksize=(50, 50, 480, 20), meta=np.ndarray>
    u_vvisc        (ocean_time, s_rho, lat_rho, lon_u) float32 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    v_vvisc        (ocean_time, s_rho, lat_v, lon_rho) float32 dask.array<chunksize=(50, 50, 480, 20), meta=np.ndarray>
    u_hvisc        (ocean_time, s_rho, lat_rho, lon_u) float32 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    v_hvisc        (ocean_time, s_rho, lat_v, lon_rho) float32 dask.array<chunksize=(50, 50, 480, 20), meta=np.ndarray>
    u_xvisc        (ocean_time, s_rho, lat_rho, lon_u) float32 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    v_xvisc        (ocean_time, s_rho, lat_v, lon_rho) float32 dask.array<chunksize=(50, 50, 480, 20), meta=np.ndarray>
    u_yvisc        (ocean_time, s_rho, lat_rho, lon_u) float32 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    v_yvisc        (ocean_time, s_rho, lat_v, lon_rho) float32 dask.array<chunksize=(50, 50, 480, 20), meta=np.ndarray>
    u_accel        (ocean_time, s_rho, lat_rho, lon_u) float32 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    v_accel        (ocean_time, s_rho, lat_v, lon_rho) float32 dask.array<chunksize=(50, 50, 480, 20), meta=np.ndarray>
    temp_hadv      (ocean_time, s_rho, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    temp_xadv      (ocean_time, s_rho, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    temp_yadv      (ocean_time, s_rho, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    temp_vadv      (ocean_time, s_rho, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    temp_hdiff     (ocean_time, s_rho, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    temp_xdiff     (ocean_time, s_rho, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    temp_ydiff     (ocean_time, s_rho, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    temp_sdiff     (ocean_time, s_rho, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    temp_vdiff     (ocean_time, s_rho, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    temp_rate      (ocean_time, s_rho, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    salt_hadv      (ocean_time, s_rho, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    salt_xadv      (ocean_time, s_rho, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    salt_yadv      (ocean_time, s_rho, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    salt_vadv      (ocean_time, s_rho, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    salt_hdiff     (ocean_time, s_rho, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    salt_xdiff     (ocean_time, s_rho, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    salt_ydiff     (ocean_time, s_rho, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    salt_sdiff     (ocean_time, s_rho, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    salt_vdiff     (ocean_time, s_rho, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    salt_rate      (ocean_time, s_rho, lat_rho, lon_rho) float32 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    pm_v           (lat_v, lon_rho) float64 dask.array<chunksize=(480, 20), meta=np.ndarray>
    pn_u           (lat_rho, lon_u) float64 dask.array<chunksize=(481, 19), meta=np.ndarray>
    pm_u           (lat_rho, lon_u) float64 dask.array<chunksize=(481, 19), meta=np.ndarray>
    pn_v           (lat_v, lon_rho) float64 dask.array<chunksize=(480, 20), meta=np.ndarray>
    pm_psi         (lat_v, lon_u) float64 dask.array<chunksize=(480, 19), meta=np.ndarray>
    pn_psi         (lat_v, lon_u) float64 dask.array<chunksize=(480, 19), meta=np.ndarray>
    dx             (lat_rho, lon_rho) float64 dask.array<chunksize=(481, 20), meta=np.ndarray>
    dx_u           (lat_rho, lon_u) float64 dask.array<chunksize=(481, 19), meta=np.ndarray>
    dx_v           (lat_v, lon_rho) float64 dask.array<chunksize=(480, 20), meta=np.ndarray>
    dx_psi         (lat_v, lon_u) float64 dask.array<chunksize=(480, 19), meta=np.ndarray>
    dy             (lat_rho, lon_rho) float64 dask.array<chunksize=(481, 20), meta=np.ndarray>
    dy_u           (lat_rho, lon_u) float64 dask.array<chunksize=(481, 19), meta=np.ndarray>
    dy_v           (lat_v, lon_rho) float64 dask.array<chunksize=(480, 20), meta=np.ndarray>
    dy_psi         (lat_v, lon_u) float64 dask.array<chunksize=(480, 19), meta=np.ndarray>
    dz             (ocean_time, s_rho, lat_rho, lon_rho) float64 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    dz_w           (ocean_time, s_w, lat_rho, lon_rho) float64 dask.array<chunksize=(50, 1, 481, 20), meta=np.ndarray>
    dz_u           (ocean_time, s_rho, lat_rho, lon_u) float64 dask.array<chunksize=(50, 50, 481, 19), meta=np.ndarray>
    dz_w_u         (ocean_time, s_w, lat_rho, lon_u) float64 dask.array<chunksize=(50, 1, 481, 19), meta=np.ndarray>
    dz_v           (ocean_time, s_rho, lat_v, lon_rho) float64 dask.array<chunksize=(50, 50, 480, 20), meta=np.ndarray>
    dz_w_v         (ocean_time, s_w, lat_v, lon_rho) float64 dask.array<chunksize=(50, 1, 480, 20), meta=np.ndarray>
    dz_psi         (ocean_time, s_rho, lat_v, lon_u) float64 dask.array<chunksize=(50, 50, 480, 19), meta=np.ndarray>
    dz_w_psi       (ocean_time, s_w, lat_v, lon_u) float64 dask.array<chunksize=(50, 1, 480, 19), meta=np.ndarray>
    dz0            (s_rho, lat_rho, lon_rho) float64 dask.array<chunksize=(50, 481, 20), meta=np.ndarray>
    dz_w0          (s_w, lat_rho, lon_rho) float64 dask.array<chunksize=(1, 481, 20), meta=np.ndarray>
    dz_u0          (s_rho, lat_rho, lon_u) float64 dask.array<chunksize=(50, 481, 19), meta=np.ndarray>
    dz_w_u0        (s_w, lat_rho, lon_u) float64 dask.array<chunksize=(1, 481, 19), meta=np.ndarray>
    dz_v0          (s_rho, lat_v, lon_rho) float64 dask.array<chunksize=(50, 480, 20), meta=np.ndarray>
    dz_w_v0        (s_w, lat_v, lon_rho) float64 dask.array<chunksize=(1, 480, 20), meta=np.ndarray>
    dz_psi0        (s_rho, lat_v, lon_u) float64 dask.array<chunksize=(50, 480, 19), meta=np.ndarray>
    dz_w_psi0      (s_w, lat_v, lon_u) float64 dask.array<chunksize=(1, 480, 19), meta=np.ndarray>
    dA             (lat_rho, lon_rho) float64 dask.array<chunksize=(481, 20), meta=np.ndarray>
    dA_u           (lat_rho, lon_u) float64 dask.array<chunksize=(481, 19), meta=np.ndarray>
    dA_v           (lat_v, lon_rho) float64 dask.array<chunksize=(480, 20), meta=np.ndarray>
    dA_psi         (lat_v, lon_u) float64 dask.array<chunksize=(480, 19), meta=np.ndarray>
    dV             (ocean_time, s_rho, lat_rho, lon_rho) float64 dask.array<chunksize=(50, 50, 481, 20), meta=np.ndarray>
    dV_w           (ocean_time, s_w, lat_rho, lon_rho) float64 dask.array<chunksize=(50, 1, 481, 20), meta=np.ndarray>
    dV_u           (ocean_time, s_rho, lat_rho, lon_u) float64 dask.array<chunksize=(50, 50, 481, 19), meta=np.ndarray>
    dV_w_u         (ocean_time, s_w, lat_rho, lon_u) float64 dask.array<chunksize=(50, 1, 481, 19), meta=np.ndarray>
    dV_v           (ocean_time, s_rho, lat_v, lon_rho) float64 dask.array<chunksize=(50, 50, 480, 20), meta=np.ndarray>
    dV_w_v         (ocean_time, s_w, lat_v, lon_rho) float64 dask.array<chunksize=(50, 1, 480, 20), meta=np.ndarray>
    dV_psi         (ocean_time, s_rho, lat_v, lon_u) float64 dask.array<chunksize=(50, 50, 480, 19), meta=np.ndarray>
    dV_w_psi       (ocean_time, s_w, lat_v, lon_u) float64 dask.array<chunksize=(50, 1, 480, 19), meta=np.ndarray>
Attributes:
    file:              /glade/scratch/dwhitt/tpos20/OUT/ocean_dia_tpos20_3_d2...
    format:            netCDF-3 64bit offset file
    Conventions:       CF-1.4, SGRID-0.3
    type:              ROMS/TOMS diagnostics file
    title:             TPOS 1/20deg
    var_info:          /glade/p/nsc/ncgd0043/tpos_ROMS_runs_setup_scripts/tpo...
    rst_file:          /glade/scratch/dwhitt/tpos20/OUT/ocean_rst_tpos20_3_d2...
    avg_file:          /glade/scratch/dwhitt/tpos20/OUT/ocean_avg_tpos20_3_d2...
    dia_file:          /glade/scratch/dwhitt/tpos20/OUT/ocean_dia_tpos20_3_d2...
    sta_file:          /glade/scratch/dwhitt/tpos20/OUT/ocean_sta_tpos20_3_d2...
    grd_file:          /glade/p/nsc/ncgd0043/tpos_ROMS_runs_setup_scripts/tpo...
    ini_file:          /glade/p/nsc/ncgd0043/tpos20/OUT/ocean_rst_tpos20_3_d2.nc
    frc_file_01:       /glade/p/nsc/ncgd0043/tpos_ROMS_runs_setup_scripts/tpo...
    frc_file_02:       /glade/p/nsc/ncgd0043/tpos_ROMS_runs_setup_scripts/tpo...
    frc_file_03:       /glade/p/nsc/ncgd0043/tpos_ROMS_runs_setup_scripts/tpo...
    bry_file:          /glade/p/nsc/ncgd0043/tpos_ROMS_runs_setup_scripts/tpo...
    script_file:       ./tpos20_3_diag_2.in
    spos_file:         /glade/p/nsc/ncgd0043/tpos_ROMS_runs_setup_scripts/tpo...
    NLM_LBC:           \nEDGE:  WEST   SOUTH  EAST   NORTH  \nzeta:  Che    C...
    svn_url:           https:://myroms.org/svn/src
    svn_rev:           Unversioned directory
    code_dir:          /glade/p/nsc/ncgd0043/tpos_ROMS_runs_setup_scripts/ROM...
    header_dir:        /glade/p/nsc/ncgd0043/tpos_ROMS_runs_setup_scripts/tpo...
    header_file:       tpos20.h
    os:                Linux
    cpu:               x86_64
    compiler_system:   ifort
    compiler_command:  /glade/u/apps/ch/opt/ncarcompilers/0.5.0/intel/18.0.5/...
    compiler_flags:    -fp-model precise -ip -O3
    tiling:            018x012
    history:           ROMS/TOMS, Version 3.7, Tuesday - July 30, 2019 -  5:5...
    ana_file:          ROMS/Functionals/ana_btflux.h
    CPP_options:       TPOS20, ALBEDO_JRA, ANA_BSFLUX, ANA_BTFLUX, ASSUMED_SH...
t0 = subsetter(ds).ocean_time[1].data  # pd.Timestamp("1985-Oct-02 00:00:00").to_numpy()
subds = ds.pipe(subsetter).pipe(interpolate, newz).rename(renamer)

RHS forcing#

rhs = xr.Dataset()

rhs_tracer = ["hadv", "vadv", "hdiff"]
rhs_mom = ["prsgrd", "hadv", "vadv", "hvisc"]

for var, terms in zip(
    ["u", "v", "temp", "salt"], [rhs_mom, rhs_mom, rhs_tracer, rhs_tracer]
):
    varnames = [f"{var}_{term}" for term in terms]
    rhs[var] = (subds[varnames].to_array("term").sum("term")) * 86400
    rhs[var].attrs["description"] = f"sum over {varnames}"
rhs
<xarray.Dataset>
Dimensions:      (ocean_time: 20, z_rho: 216)
Coordinates:
    nDIA         int32 180
    LwSrc        int32 0
    Zos          float64 0.02
    Fgamma       float64 0.284
    lon_u        float64 -140.0
    Zob          float64 0.02
    ntsDIA       int32 1
    rdrg         float64 0.000265
    mask_rho     float64 1.0
    z_rho_psi    (ocean_time, z_rho) float64 -107.8 -107.3 ... -0.7773 -0.2768
    gamma2       float64 1.0
    f            float64 4.514e-18
    ndefAVG      int32 0
    dt           float64 120.0
    z_rho_v      (ocean_time, z_rho) float64 -107.7 -107.2 ... -0.776 -0.2764
    ntimes       int32 25201
    lon_rho      float64 -140.0
    ndtfast      int32 20
    mask_u       float64 1.0
    lon_v        float64 -140.0
    Falpha       float64 2.0
    M2nudg       float64 0.0
  * z_rho        (z_rho) float64 -107.8 -107.2 -106.8 ... -1.25 -0.75 -0.25
    LuvSrc       int32 0
    dstart       object 0027-09-02 00:00:00
    ntsAVG       int32 1
    xl           float64 8.163e+06
    Vtransform   int32 1
    dtfast       float64 6.0
    Lm2PSH       int32 0
    pn           float64 0.0001799
    LsshCLM      int32 0
    LnudgeM2CLM  int32 0
    mask_v       float64 1.0
    rdrg2        float64 0.003
    z_rho0       (z_rho) float64 -108.1 -107.6 -107.1 ... -1.577 -1.077 -0.5769
    lat_u        float64 1.772e-12
    LnudgeM3CLM  int32 0
    Cs_r         (z_rho) float64 -0.02114 -0.02104 ... -0.0002008 -0.0001072
    rho0         float64 1.025e+03
    hc           float64 75.0
    nHIS         int32 3600
    ndefDIA      int32 0
    el           float64 2.674e+06
    angle        float64 4.106e-14
    pm           float64 0.0001799
    Znudg        float64 0.0
    spherical    int32 1
    ndefHIS      int32 0
    Akv_bak      float64 1e-05
    z_rho_v0     (z_rho) float64 -108.0 -107.5 -107.0 ... -1.576 -1.076 -0.5765
    nl_visc4     float64 1.24e+09
    grid         int32 1
    z_rho_u0     (z_rho) float64 -108.2 -107.7 -107.2 ... -1.578 -1.078 -0.5774
    nSTA         int32 30
    M3nudg       float64 1.0
    nRST         int32 3600
    Tcline       float64 75.0
    Lm2CLM       int32 0
    Lm3PSH       int32 0
    mask_psi     float64 1.0
    Lm3CLM       int32 0
    lat_psi      float64 -0.025
    s_rho        (z_rho) float64 -0.2463 -0.2455 -0.2446 ... -0.003013 -0.001633
    z_rho_u      (ocean_time, z_rho) float64 -107.8 -107.3 ... -0.777 -0.2766
  * ocean_time   (ocean_time) datetime64[ns] 1985-10-02T03:00:00 ... 1985-10-...
    lat_rho      float64 1.45e-14
    LuvSponge    int32 0
    nAVG         int32 180
    lat_v        float64 -0.025
    h            float64 4.314e+03
    z_rho_psi0   (z_rho) float64 -108.2 -107.7 -107.2 ... -1.578 -1.078 -0.5774
    lon_psi      float64 -140.0
    Vstretching  int32 1
Data variables:
    u            (ocean_time, z_rho) float64 0.0935 0.09244 ... 0.1627 0.1625
    v            (ocean_time, z_rho) float64 0.106 0.1083 ... -0.0536 -0.05361
    temp         (ocean_time, z_rho) float64 -0.8086 -0.7765 ... -0.02899
    salt         (ocean_time, z_rho) float64 -0.02245 -0.02013 ... 0.0006571
write_to_txt(
    rhs.rename(renamer).pipe(interleave_time, t0),
    outdir=outdir,
    prefix="ocean_colpsh_",
)

initial condition#

substate = (subds[["U", "W", "T", "S"]].compute()) - offset

ic = substate.cf.sel(time=t0)
write_to_txt(ic.compute(), outdir=outdir, prefix="ocean_ic_")

restoring#

write_to_txt(interleave_time(substate, t0), outdir=outdir, prefix="ocean_colfrc_")

surface forcing#

frc = subds[["sustr", "svstr", "swrad", "ssflux"]]
frc["shfluxminusswrad"] = subds["shflux"] - frc["swrad"]

time = subds.indexes["ocean_time"].to_numpy()
frc["time"] = ("ocean_time", (time - t0).astype("timedelta64[s]").astype(float))
frc.load()

write_to_txt(frc, outdir=outdir, prefix="ocean_frc_")
# Test flux calculations
shflux2 = subds[["lwrad", "latent", "sensible"]].to_array().sum("variable")

subds.swrad.plot()
subds.shflux.plot()
(-subds.swrad + subds.shflux).plot(marker="x")
(shflux2).plot()
[<matplotlib.lines.Line2D at 0x2b245c327940>]
_images/a0389fd513ea038aefe1b54a593931c3b2753f19667d818efb821ee01021464d.png