Tracer budget closure in MOM6

Contributors: Graeme MacGilchrist, Stephen Griffies

Notes: Much of the nomenclature here is taken from the MOM5 documentation, which can be cited as Griffies et al. (2012) Elements of the Modular Ocean Model, NOAA/GFDL Technical Report.

In this tutorial, we outline the terms required to close the heat, salt and thickness budgets in MOM6. We show this for diagnostics that are output both on the native grid and remapped to a diagnostic grid (\(z^*\) and \(\rho_2\)).

Conservation of vertically-extensive tracer content

Conservation of tracer concentration \(\phi\) can be written

\[\rho\frac{D\phi}{Dt} = -\nabla\cdot \mathbf{J}_\phi + \rho S_{\phi}\]

where \(D/Dt\) is the material derivative, \(\mathbf{J}_\phi\) is the diffusive (non-advective) tracer flux, and \(S_{\phi}\) are sources/sinks.

Due to its Lagrangian vertical coordinate, MOM6 is formulated to conserve the vertically-extensive (concentration times thickness) tracer content in each layer at each point on the horizontal grid. Integrated vertically in a model layer, the material derivative in the above expression takes the semi-continuous form:

\[\int^{z_{k-1}}_{z_k}\Big(\rho\frac{D\phi}{Dt}\Big)\,dz = \partial_t (\rho \phi h) + \nabla_s \cdot (\rho \phi \mathbf{u}\, h) + \Delta_k(\rho \phi w^{(s)})\,,\]

where \(\partial_t\) is the grid-cell time derivative, \(h\) is the layer thickness, \(\nabla_s\) is the gradient operator parallel to the general vertical coordinate \(s\), \(\mathbf{u}\) is the horizontal velocity within a layer, \(\Delta_k\) is a discrete difference operator between the upper and lower bounds of the layer, and \(w^{(s)}\) is the across-layer velocity. Setting the tracer concentration, \(\phi\), to 1, this becomes an expression for mass conservation (equivalent to thickness conservation under the Boussinesq approximation and with fixed horizontal grid spacing).

In closing the MOM6 budget, the first term on the right hand side of this expression corresponds to the total tendency of the tracer content per horizontal area in a grid cell, the second term corresponds to the tendency due to the divergence of along-layer advection, and the third term corresponds to the tendency due to the vertical remapping.

Correspondingly, the vertically-integrated right hand side of the tracer concentration conservation equation can be written

\[\int^{z_{k-1}}_{z_k}(-\nabla\cdot \mathbf{J}_\phi)\,dz = -\nabla_s\cdot(\mathbf{J}_\phi\,h)-\Delta_k(\mathbf{J}^{(s)}_\phi) + \rho S_{\phi}\,h,\]

where \(\mathbf{J}^{(s)}_\phi\) is the vertical component of the tracer flux.

Non-local isoneutral diffusion and the challenge of closing budgets from fluxes

Isoneutral diffusion in the arbitrary vertical coordinate of MOM6 is handled in such a way that means tracer fluxes can be “non-local”. That is, a diffusive flux of tracer for example in the \(x\)-direction out of layer \(k\) at horizontal grid-point \((i,j)\) need not go directly into layer \(k\) at horizontal grid-point \((i+1,j)\). Rather, it can go into any layer at grid-point \((i+1,j)\) that satifies the conditions for iso-neutrality.

This makes it, at present, impossible to fully close layer-wise tracer budgets using the diffusive fluxes at a grid-cell’s lateral boundaries, since they are closed only in a column-integrated sense. Instead, we can only fully close the layerwise budgets from the tendencies of the tracer content for a grid cell. That is because, even the keeping track of lateral isoneutral fluxes is impossible, we can still calculate the total tendency due to these fluxes at each point and in each layer.

Note the that colum-integrated tracer budget can be closed from the lateral flux and source terms directly.

Terms in tracer budget

The processes influencing the heat, salt, and thickness budgets in MOM6, and the standard diagnostic names for their corresponding tendencies are given in the following table:

Process	Heat	Salt	Thickness
Grid-cell time tendency	opottemptend	osalttend	dhdt
Within-layer horizontal advection	T_advection_xy	S_advection_xy	dynamics_h_tendency
Across-layer advection (remapping)	Th_tendency_vert_remap	Sh_tendency_vert_remap	vert_remap_h_tendency
Boundary forcing	boundary_forcing_heat_tendency	boundary_forcing_salt_tendency	boundary_forcing_h_tendency
Neutral diffusion	opottemppmdiff	osaltpmdiff	N/A
Across-layer diffusion	opottempdiff	osaltdiff	N/A
Internal (geothermal) heating	internal_heat_heat_tendency	N/A	N/A
Frazil ice formation	frazil_heat_tendency	N/A	N/A

Further separation of the budget terms (e.g. separating boundary forcing of heat into components due to shortwave and longwave radiation, latent and sensible heat fluxes) is possible. Indeed many of these diagnostics already exist and can be found in any available_diags file. Here we are presenting just the most straightforward budget closure.

Sign convention

The sign convention is such that the grid-cell time tendency is equal to the sum of all the other terms. As such, in equating the left and right sides of the tracer conservation equation, it is necessary to take away the along-layer and across-layer advective terms from the grid-cell time tendency.

So, for example, dhdt - dynamics_h_tendency - vert_remap_h_tendency = boundary_forcing_h_tendency

Showing budget closure

Here, we confirm that, given the terms in the table above, the left hand side (material derivative) and right hand side (diffusive fluxes and source/sink) of the tracer conservation equation are equal and thus sum to zero (to machine precision). We confirm this for (a) the full depth integral, (b) a single profile, and (c) integrated within each layer (full 3D closure).

For this purpose we look at a simple, Baltic Sea configuration of OM4, available here. Details on downloading and compiling the code, and runnig the example simulations can be found on the MOM6-examples wiki. We ran the simulation for one year, outputting diagnostics as monthly averages, and here just look at the first month.

# Import packages
import xarray as xr
from matplotlib import pyplot as plt
# Don't display filter warnings
import warnings
warnings.filterwarnings("ignore")
# Set figure font size
plt.rcParams.update({'font.size':14})

# Specify budget terms
terms = {}
terms['heat'] = ['opottemptend','T_advection_xy','Th_tendency_vert_remap',
                 'boundary_forcing_heat_tendency','opottempdiff','opottemppmdiff',
                 'frazil_heat_tendency','internal_heat_heat_tendency']
terms['salt'] = ['osalttend','S_advection_xy','Sh_tendency_vert_remap',
                 'boundary_forcing_salt_tendency','osaltdiff','osaltpmdiff']
terms['thickness'] = ['dhdt','dynamics_h_tendency','vert_remap_h_tendency',
                      'boundary_forcing_h_tendency']
# Specify different diagnostic grids
grids = ['native','z','rho2']
# Specify vertical index name in grid
vertind = {'native':'zl','z':'z_l','rho2':'rho2_l'}

# Define a function to sum up the terms on the LHS and RHS of the tracer
# conservation equation
def calc_budget(ds,terms):
    lhs = xr.Dataset()
    rhs = xr.Dataset()
    for key in terms.keys():
        lhs[key] = xr.zeros_like(ds[terms[key][0]])
        for term in terms[key][0:3]:
            # Flip the sign for the advective tendencies
            if term in ['opottemptend','osalttend','dhdt']:
                sign=1
            else:
                sign=-1
            lhs[key] += sign*ds[term]
        rhs[key] = xr.zeros_like(ds[terms[key][3]])
        for term in terms[key][3:]:
            sign=1
            rhs[key] += sign*ds[term]
    return lhs, rhs

# Load data on the native grid, z* grid and rho2 grids
# There are separate files for each grid,
# so loop through to load and place each dataset into a dictionary
# Within the loop, calculate the LHS and RHS for each tracer for each grid
rootdir = '/archive/gam/MOM6-examples/ice_ocean_SIS2/Baltic_OM4_025/1yr/'
prefix = '19000101.ocean_'
ds = {}
lhs = {}
rhs = {}
for grid in grids:
    # Diagnostics were saved into different files
    suffixs = ['thck','heat','salt']
    ds[grid] = xr.Dataset()
    for suffix in suffixs:
        filename = prefix+grid+'_'+suffix+'*.nc'
        dsnow = xr.open_mfdataset(rootdir+filename).isel(time=0)
        ds[grid] = xr.merge([ds[grid],dsnow])
    # Calculate the two sides of the tracer equation on the native grid
    lhs[grid],rhs[grid] = calc_budget(ds[grid],terms)

Diagnostics output on the native grid

MOM6 has the capability of outputting diagnostics either on the native grid (the grid on which the equations are solved) or a user-defined diagnostic grid, where the diagnostics have been remapped onto, for example, \(z^*\) or \(\rho_2\) levels.

Here we show the budget closure on the native grid.

Full-depth integral budget closure

Below, we plot the vertical integral of the sides of the tracer conservation equation, and their difference.

Variables were saved at double precision, so for the budget to be considered “closed” the difference between the LHS and RHS should be approximately 12 orders of magnitude smaller than the contributing terms. As can be seen below, this is true for each of heat, salt and thickness.

grid = 'native'
# Plot the vertical sum
fig,ax = plt.subplots(figsize=(16,8),nrows=3,ncols=3)
count=0
for key in terms.keys():
    # LHS
    im = ax[0,count].pcolormesh(
        lhs[grid][key].sum(vertind[grid]).squeeze())
    plt.colorbar(im,ax=ax[0,count])
    ax[0,count].set_title(key)
    # RHS
    im = ax[1,count].pcolormesh(
        rhs[grid][key].sum(vertind[grid]).squeeze())
    plt.colorbar(im,ax=ax[1,count])
    # Difference
    im = ax[2,count].pcolormesh(
        (lhs[grid][key]-rhs[grid][key]).sum(vertind[grid]).squeeze(),
        vmin=0)
    plt.colorbar(im,ax=ax[2,count])
    count+=1
ax[0,0].set_ylabel('LHS')
ax[1,0].set_ylabel('RHS')
ax[2,0].set_ylabel('difference')
plt.tight_layout()

Local (cell-by-cell) budget closure

It is important to check that the local (cell-by-cell) budget closes, not just that of the full-depth integral. This is because certain processes in the MOM6 algorithm, such as vertical remapping and vertical diffusion, sum to zero in the vertical by construction. Thus, it is possible that those terms could be diagnosed incorrectly, while the full-depth budget remained closed.

We can check that the budget is closed at individual grid cells by looking at a profile at a single location. In this case, it makes sense to define the y-axis as the cumulative thickness of the layers, which defines an approximate depth, and allows like-for-like comparison (particularly when plotting output from different diagnostic grids, as below).

As before, the difference between the sides of the tracer conservation equation is extremely small compared with either side.

# Plot a profile
# [In this case, plot the terms against the vertical sum of the cell
# thicknesses, so that the y-axis is depth]
lon=5
lat=56
select = {'xh':lon,'yh':lat}
fig,ax = plt.subplots(figsize=(9,6),ncols=3)
count=0
z = ds[grid]['thkcello'].sel(select,method='nearest').squeeze().cumsum(vertind[grid])
for key in terms.keys():
    ax[count].plot(
        lhs[grid][key].sel(select,method='nearest').squeeze(),z,label='LHS')
    ax[count].plot(
        rhs[grid][key].sel(select,method='nearest').squeeze(),z,label='RHS')
    ax[count].plot(
        (lhs[grid][key]-rhs[grid][key]).sel(select,method='nearest').squeeze(),z,label='difference')
    ax[count].invert_yaxis()
    ax[count].set_title(key)
    count+=1
fig.suptitle('lon: '+str(lon)+', lat: '+str(lat))
ax[2].legend()

<matplotlib.legend.Legend at 0x2ae9607febe0>

Root-mean-square error

In fact, the most comprehensive way to confirm that the budget is closed is to examine the root-mean-square error (difference) between the sides of the tracer conservation equation. This could be done in any of the \(x\), \(y\), or \(z\), but here we do it in the vertical:

\[RMSE(i,j) = \sqrt{\sum^{n_k}_{k=1}(LHS(i,j,k)-RHS(i,j,k))^2/n_k}\,,\]

where \(k\) is the vertical index.

Again, the RMSE should be around 12 orders of magnitude less than either the left or right hand side terms. By comparison to the vertical sum figure above, it is clear that this is the case.

# Calculate the RMSE
rmse = xr.ufuncs.sqrt(xr.ufuncs.square(lhs[grid]-rhs[grid]).mean('zl')).squeeze()
fig,ax = plt.subplots(figsize=(16,4),ncols=3)
count=0
for key in terms.keys():
    im = ax[count].pcolormesh(rmse[key])
    plt.colorbar(im,ax=ax[count])
    ax[count].set_title(key)
    count+=1
plt.tight_layout()

Diagnostics output on a user-defined vertical grid

Outputting variables on a diagnostic grid is done through a remapping procedure, and it is important to confirm that the procedure does not break the budget closure.

[TO ADD: More detail on the remapping to diagnostic grids]

Local (cell-by-cell) budget closure

We show here a direct, cell-by-cell comparison between the difference arising from diagnostics on the native grid and on the user-defined grids. As above, we plot variables against the cumulative thickness, such that they are defined at the same depth, irrespective of the diagnostic grid.

Here, we can see that the errors are all fairly similary between the different grids and all very small (by comparison to the terms in the profile figure above), implying that the budgets are equivalently closed on the native and diagnostic grids.

lon=5
lat=90
select = {'xh':lon,'yh':lat}

fig,ax = plt.subplots(figsize=(9,6),ncols=3)
for grid in grids:
    z = ds[grid]['thkcello'].sel(select,method='nearest').squeeze().cumsum(vertind[grid])
    count=0
    for key in terms.keys():
        ax[count].plot(
            (lhs[grid][key]-rhs[grid][key]).sel(select,method='nearest').squeeze(),z,label=grid)
        ax[count].set_title(key)
        count+=1
ax[0].invert_yaxis()
ax[1].invert_yaxis()
ax[2].invert_yaxis()
ax[2].legend()

<matplotlib.legend.Legend at 0x2ae96c21c9a0>

Root-mean-square error

Finally, we compare the RMSE for each of the output grids. On each of the diagnostic grids, the RMSE is small, and comparable to that on the native grid. The errors are very marginally greater on the \(z^*\) grid, but still a tiny fraction of the contributing terms.

# Plot the vertical sum
fig,ax = plt.subplots(figsize=(16,8),nrows=3,ncols=3)
count_grid=0
for grid in grids:
    # Calculate RMSE
    rmse = xr.ufuncs.sqrt(
        xr.ufuncs.square(lhs[grid]-rhs[grid]).mean(vertind[grid])).squeeze()
    count_key=0
    for key in terms.keys():
        # RMSE
        im = ax[count_grid,count_key].pcolormesh(rmse[key].squeeze())
        plt.colorbar(im,ax=ax[count_grid,count_key])
        ax[count_grid,count_key].set_title(grid+', '+key)
        count_key+=1
    count_grid+=1
plt.tight_layout()

Example diag_table

The diag_table used for this Baltic Sea configuration, including additional diagnostics necessary for the calculation of watermass transformation (see Watermass transformation in MOM6), is shown below.

Baltic_OM4_025
1900 1 1 0 0 0

#====================================================================================================
# MOM6
#====================================================================================================
# MOM6 ocean diagnostics files

"ocean_native_thck%4yr_%2mo",     1, "months",   1, "days", "time", 1, "months"
"ocean_native_heat%4yr_%2mo",     1, "months",   1, "days", "time", 1, "months"
"ocean_native_salt%4yr_%2mo",     1, "months",   1, "days", "time", 1, "months"
"ocean_native_snap%4yr_%2mo",     1, "months",   1, "days", "time", 1, "months"
"ocean_rho2_thck%4yr_%2mo",       1, "months",   1, "days", "time", 1, "months"
"ocean_rho2_heat%4yr_%2mo",       1, "months",   1, "days", "time", 1, "months"
"ocean_rho2_salt%4yr_%2mo",       1, "months",   1, "days", "time", 1, "months"
"ocean_rho2_snap%4yr_%2mo",       1, "months",   1, "days", "time", 1, "months"
"ocean_z_thck%4yr_%2mo",          1, "months",   1, "days", "time", 1, "months"
"ocean_z_heat%4yr_%2mo",          1, "months",   1, "days", "time", 1, "months"
"ocean_z_salt%4yr_%2mo",          1, "months",   1, "days", "time", 1, "months"
"ocean_z_snap%4yr_%2mo",          1, "months",   1, "days", "time", 1, "months"
"ocean_static",          -1, "hours",   1, "days", "time" # ocean_static is a protected name. Do not change this line.

# --------------------------- #
# Static grid variables       #
# --------------------------- #
"ocean_model", "areacello",   "areacello",   "ocean_static", "all", "none", "none", 2
"ocean_model", "geolon",      "geolon",      "ocean_static", "all", "none", "none", 2
"ocean_model", "geolat",      "geolat",      "ocean_static", "all", "none", "none", 2
"ocean_model", "geolon_c",    "geolon_c",    "ocean_static", "all", "none", "none", 2
"ocean_model", "geolat_c",    "geolat_c",    "ocean_static", "all", "none", "none", 2
"ocean_model", "geolon_u",    "geolon_u",    "ocean_static", "all", "none", "none", 2
"ocean_model", "geolat_u",    "geolat_u",    "ocean_static", "all", "none", "none", 2
"ocean_model", "geolon_v",    "geolon_v",    "ocean_static", "all", "none", "none", 2
"ocean_model", "geolat_v",    "geolat_v",    "ocean_static", "all", "none", "none", 2
"ocean_model", "area_t",      "area_t",      "ocean_static", "all", "none", "none", 2
"ocean_model", "area_u",      "area_u",      "ocean_static", "all", "none", "none", 2
"ocean_model", "area_v",      "area_v",      "ocean_static", "all", "none", "none", 2
"ocean_model", "area_q",      "area_q",      "ocean_static", "all", "none", "none", 2
"ocean_model", "depth_ocean", "depth_ocean", "ocean_static", "all", "none", "none", 2
"ocean_model", "wet",         "wet",         "ocean_static", "all", "none", "none", 2
"ocean_model", "wet_c",       "wet_c",       "ocean_static", "all", "none", "none", 2
"ocean_model", "wet_u",       "wet_u",       "ocean_static", "all", "none", "none", 2
"ocean_model", "wet_v",       "wet_v",       "ocean_static", "all", "none", "none", 2
"ocean_model", "Coriolis",    "Coriolis",    "ocean_static", "all", "none", "none", 2
"ocean_model", "dxt",         "dxt",         "ocean_static", "all", "none", "none", 2
"ocean_model", "dyt",         "dyt",         "ocean_static", "all", "none", "none", 2
"ocean_model", "dxCu",        "dxCu",        "ocean_static", "all", "none", "none", 2
"ocean_model", "dyCu",        "dyCu",        "ocean_static", "all", "none", "none", 2
"ocean_model", "dxCv",        "dxCv",        "ocean_static", "all", "none", "none", 2
"ocean_model", "dyCv",        "dyCv",        "ocean_static", "all", "none", "none", 2
"ocean_model", "dxCvo",       "dxCvo",       "ocean_static", "all", "none", "none", 2
"ocean_model", "dyCuo",       "dyCuo",       "ocean_static", "all", "none", "none", 2

##############
### NATIVE ###
# ------------------------------------------------- #
# Heat + budget terms [get from Andrew Shao github] #
# ------------------------------------------------- #

# MEAN
"ocean_model",  "temp",                         "temp"                                  "ocean_native_heat%4yr_%2mo", "all", "mean", "none", 2

"ocean_model",  "opottemptend",                 "opottemptend",                         "ocean_native_heat%4yr_%2mo", "all", "mean", "none", 1
"ocean_model",  "opottempdiff",                 "opottempdiff",                         "ocean_native_heat%4yr_%2mo", "all", "mean", "none", 1
"ocean_model",  "opottemppmdiff",               "opottemppmdiff",                       "ocean_native_heat%4yr_%2mo", "all", "mean", "none", 1
"ocean_model",  "boundary_forcing_heat_tendency","boundary_forcing_heat_tendency",      "ocean_native_heat%4yr_%2mo", "all", "mean", "none", 1
"ocean_model",  "frazil_heat_tendency",         "frazil_heat_tendency",                 "ocean_native_heat%4yr_%2mo", "all", "mean", "none", 1
"ocean_model",  "T_advection_xy",               "T_advection_xy",                       "ocean_native_heat%4yr_%2mo", "all", "mean", "none", 1
"ocean_model",  "Th_tendency_vert_remap",       "Th_tendency_vert_remap",               "ocean_native_heat%4yr_%2mo", "all", "mean", "none", 1
"ocean_model",  "internal_heat_heat_tendency",  "internal_heat_heat_tendency",          "ocean_native_heat%4yr_%2mo", "all", "mean", "none", 1

# SNAP
"ocean_model",  "temp",                         "temp"                                  "ocean_native_snap%4yr_%2mo", "all", "none", "none", 2

# ------------------- #
# Salt + budget terms #
# ------------------- #

"ocean_model",  "salt",                         "salt",                                 "ocean_native_salt%4yr_%2mo", "all", "mean", "none", 2

"ocean_model",  "osalttend",                    "osalttend",                            "ocean_native_salt%4yr_%2mo", "all", "mean", "none", 1
"ocean_model",  "osaltdiff",                    "osaltdiff",                            "ocean_native_salt%4yr_%2mo", "all", "mean", "none", 1
"ocean_model",  "osaltpmdiff",                  "osaltpmdiff",                          "ocean_native_salt%4yr_%2mo", "all", "mean", "none", 1
"ocean_model",  "boundary_forcing_salt_tendency","boundary_forcing_salt_tendency",      "ocean_native_salt%4yr_%2mo", "all", "mean", "none", 1
"ocean_model",  "S_advection_xy",               "S_advection_xy",                       "ocean_native_salt%4yr_%2mo", "all", "mean", "none", 1
"ocean_model",  "Sh_tendency_vert_remap",       "Sh_tendency_vert_remap",               "ocean_native_salt%4yr_%2mo", "all", "mean", "none", 1

# SNAP
"ocean_model",  "salt",                         "salt"                                  "ocean_native_snap%4yr_%2mo", "all", "none", "none", 2

# ------------------- #
# Thickness + budget terms #
# ------------------- #

# MEAN
"ocean_model",  "thkcello",                     "thkcello",                             "ocean_native_thck%4yr_%2mo", "all", "mean", "none", 2

"ocean_model",  "dhdt",                         "dhdt",                                 "ocean_native_thck%4yr_%2mo", "all", "mean", "none", 1
"ocean_model",  "dynamics_h_tendency",          "dynamics_h_tendency",                  "ocean_native_thck%4yr_%2mo", "all", "mean", "none", 1
"ocean_model",  "vert_remap_h_tendency",        "vert_remap_h_tendency",                "ocean_native_thck%4yr_%2mo", "all", "mean", "none", 1
"ocean_model",  "boundary_forcing_h_tendency",  "boundary_forcing_h_tendency",          "ocean_native_thck%4yr_%2mo", "all", "mean", "none", 1
"ocean_model",  "internal_heat_h_tendency",     "internal_heat_h_tendency",             "ocean_native_thck%4yr_%2mo", "all", "mean", "none", 1

# SNAP
"ocean_model",  "thkcello",                     "thkcello"                              "ocean_native_snap%4yr_%2mo", "all", "none", "none", 2
"ocean_model",  "e",                            "e"                                     "ocean_native_snap%4yr_%2mo", "all", "none", "none", 2

# ----------------------------------------- #
# Other terms necessary for WMT calculation #
# ----------------------------------------- #

# MEAN
"ocean_model",  "rhoinsitu",                    "rhoinsitu",                            "ocean_native_thck%4yr_%2mo", "all", "mean", "none", 2
"ocean_model",  "rhopot0",                      "rhopot0",                              "ocean_native_thck%4yr_%2mo", "all", "mean", "none", 2
"ocean_model",  "rhopot2",                      "rhopot2",                              "ocean_native_thck%4yr_%2mo", "all", "mean", "none", 2
"ocean_model",  "drho_dT",                      "drho_dT",                              "ocean_native_thck%4yr_%2mo", "all", "mean", "none", 2
"ocean_model",  "drho_dS",                      "drho_dS",                              "ocean_native_thck%4yr_%2mo", "all", "mean", "none", 2
"ocean_model",  "umo",                          "umo",                                  "ocean_native_thck%4yr_%2mo", "all", "mean", "none", 2
"ocean_model",  "vmo",                          "vmo",                                  "ocean_native_thck%4yr_%2mo", "all", "mean", "none", 2
"ocean_model",  "e",                            "e",                                    "ocean_native_thck%4yr_%2mo", "all", "mean", "none", 2

# SNAP
"ocean_model",  "rhopot2",                      "rhopot2"                               "ocean_native_snap%4yr_%2mo", "all", "none", "none", 2

############
### RHO2 ###
# ------------------------------------------------- #
# Heat + budget terms [get from Andrew Shao github] #
# ------------------------------------------------- #

# MEAN
"ocean_model_rho2",  "temp",                         "temp"                                  "ocean_rho2_heat%4yr_%2mo", "all", "mean", "none", 2

"ocean_model_rho2",  "opottemptend",                 "opottemptend",                         "ocean_rho2_heat%4yr_%2mo", "all", "mean", "none", 1
"ocean_model_rho2",  "opottempdiff",                 "opottempdiff",                         "ocean_rho2_heat%4yr_%2mo", "all", "mean", "none", 1
"ocean_model_rho2",  "opottemppmdiff",               "opottemppmdiff",                       "ocean_rho2_heat%4yr_%2mo", "all", "mean", "none", 1
"ocean_model_rho2",  "boundary_forcing_heat_tendency","boundary_forcing_heat_tendency",      "ocean_rho2_heat%4yr_%2mo", "all", "mean", "none", 1
"ocean_model_rho2",  "frazil_heat_tendency",         "frazil_heat_tendency",                 "ocean_rho2_heat%4yr_%2mo", "all", "mean", "none", 1
"ocean_model_rho2",  "T_advection_xy",               "T_advection_xy",                       "ocean_rho2_heat%4yr_%2mo", "all", "mean", "none", 1
"ocean_model_rho2",  "Th_tendency_vert_remap",       "Th_tendency_vert_remap",               "ocean_rho2_heat%4yr_%2mo", "all", "mean", "none", 1
"ocean_model_rho2",  "internal_heat_heat_tendency",  "internal_heat_heat_tendency",          "ocean_rho2_heat%4yr_%2mo", "all", "mean", "none", 1

# SNAP
"ocean_model_rho2",  "temp",                         "temp"                                  "ocean_rho2_snap%4yr_%2mo", "all", "none", "none", 2

# ------------------- #
# Salt + budget terms #
# ------------------- #

"ocean_model_rho2",  "salt",                         "salt",                                 "ocean_rho2_salt%4yr_%2mo", "all", "mean", "none", 2

"ocean_model_rho2",  "osalttend",                    "osalttend",                            "ocean_rho2_salt%4yr_%2mo", "all", "mean", "none", 1
"ocean_model_rho2",  "osaltdiff",                    "osaltdiff",                            "ocean_rho2_salt%4yr_%2mo", "all", "mean", "none", 1
"ocean_model_rho2",  "osaltpmdiff",                  "osaltpmdiff",                          "ocean_rho2_salt%4yr_%2mo", "all", "mean", "none", 1
"ocean_model_rho2",  "boundary_forcing_salt_tendency","boundary_forcing_salt_tendency",      "ocean_rho2_salt%4yr_%2mo", "all", "mean", "none", 1
"ocean_model_rho2",  "S_advection_xy",               "S_advection_xy",                       "ocean_rho2_salt%4yr_%2mo", "all", "mean", "none", 1
"ocean_model_rho2",  "Sh_tendency_vert_remap",       "Sh_tendency_vert_remap",               "ocean_rho2_salt%4yr_%2mo", "all", "mean", "none", 1

# SNAP
"ocean_model_rho2",  "salt",                         "salt"                                  "ocean_rho2_snap%4yr_%2mo", "all", "none", "none", 2

# ------------------- #
# Thickness + budget terms #
# ------------------- #

# MEAN
"ocean_model_rho2",  "thkcello",                     "thkcello",                             "ocean_rho2_thck%4yr_%2mo", "all", "mean", "none", 2

"ocean_model_rho2",  "dhdt",                         "dhdt",                                 "ocean_rho2_thck%4yr_%2mo", "all", "mean", "none", 1
"ocean_model_rho2",  "dynamics_h_tendency",          "dynamics_h_tendency",                  "ocean_rho2_thck%4yr_%2mo", "all", "mean", "none", 1
"ocean_model_rho2",  "vert_remap_h_tendency",        "vert_remap_h_tendency",                "ocean_rho2_thck%4yr_%2mo", "all", "mean", "none", 1
"ocean_model_rho2",  "boundary_forcing_h_tendency",  "boundary_forcing_h_tendency",          "ocean_rho2_thck%4yr_%2mo", "all", "mean", "none", 1
"ocean_model_rho2",  "internal_heat_h_tendency",     "internal_heat_h_tendency",             "ocean_rho2_thck%4yr_%2mo", "all", "mean", "none", 1

# SNAP
"ocean_model_rho2",  "thkcello",                     "thkcello"                              "ocean_rho2_snap%4yr_%2mo", "all", "none", "none", 2
"ocean_model_rho2",  "e",                            "e"                                     "ocean_rho2_snap%4yr_%2mo", "all", "none", "none", 2

# ----------------------------------------- #
# Other terms necessary for WMT calculation #
# ----------------------------------------- #

# MEAN
"ocean_model_rho2",  "rhoinsitu",                    "rhoinsitu",                            "ocean_rho2_thck%4yr_%2mo", "all", "mean", "none", 2
"ocean_model_rho2",  "rhopot0",                      "rhopot0",                              "ocean_rho2_thck%4yr_%2mo", "all", "mean", "none", 2
"ocean_model_rho2",  "rhopot2",                      "rhopot2",                              "ocean_rho2_thck%4yr_%2mo", "all", "mean", "none", 2
"ocean_model_rho2",  "drho_dT",                      "drho_dT",                              "ocean_rho2_thck%4yr_%2mo", "all", "mean", "none", 2
"ocean_model_rho2",  "drho_dS",                      "drho_dS",                              "ocean_rho2_thck%4yr_%2mo", "all", "mean", "none", 2
"ocean_model_rho2",  "umo",                          "umo",                                  "ocean_rho2_thck%4yr_%2mo", "all", "mean", "none", 2
"ocean_model_rho2",  "vmo",                          "vmo",                                  "ocean_rho2_thck%4yr_%2mo", "all", "mean", "none", 2
"ocean_model_rho2",  "e",                            "e",                                    "ocean_rho2_thck%4yr_%2mo", "all", "mean", "none", 2

# SNAP
"ocean_model_rho2",  "rhopot2",                      "rhopot2"                               "ocean_rho2_snap%4yr_%2mo", "all", "none", "none", 2

############
#### Z* ####
# ------------------------------------------------- #
# Heat + budget terms [get from Andrew Shao github] #
# ------------------------------------------------- #

# MEAN
"ocean_model_z",  "temp",                         "temp"                                  "ocean_z_heat%4yr_%2mo", "all", "mean", "none", 2

"ocean_model_z",  "opottemptend",                 "opottemptend",                         "ocean_z_heat%4yr_%2mo", "all", "mean", "none", 1
"ocean_model_z",  "opottempdiff",                 "opottempdiff",                         "ocean_z_heat%4yr_%2mo", "all", "mean", "none", 1
"ocean_model_z",  "opottemppmdiff",               "opottemppmdiff",                       "ocean_z_heat%4yr_%2mo", "all", "mean", "none", 1
"ocean_model_z",  "boundary_forcing_heat_tendency","boundary_forcing_heat_tendency",      "ocean_z_heat%4yr_%2mo", "all", "mean", "none", 1
"ocean_model_z",  "frazil_heat_tendency",         "frazil_heat_tendency",                 "ocean_z_heat%4yr_%2mo", "all", "mean", "none", 1
"ocean_model_z",  "T_advection_xy",               "T_advection_xy",                       "ocean_z_heat%4yr_%2mo", "all", "mean", "none", 1
"ocean_model_z",  "Th_tendency_vert_remap",       "Th_tendency_vert_remap",               "ocean_z_heat%4yr_%2mo", "all", "mean", "none", 1
"ocean_model_z",  "internal_heat_heat_tendency",  "internal_heat_heat_tendency",          "ocean_z_heat%4yr_%2mo", "all", "mean", "none", 1

# SNAP
"ocean_model_z",  "temp",                         "temp"                                  "ocean_z_snap%4yr_%2mo", "all", "none", "none", 2

# ------------------- #
# Salt + budget terms #
# ------------------- #

"ocean_model_z",  "salt",                         "salt",                                 "ocean_z_salt%4yr_%2mo", "all", "mean", "none", 2

"ocean_model_z",  "osalttend",                    "osalttend",                            "ocean_z_salt%4yr_%2mo", "all", "mean", "none", 1
"ocean_model_z",  "osaltdiff",                    "osaltdiff",                            "ocean_z_salt%4yr_%2mo", "all", "mean", "none", 1
"ocean_model_z",  "osaltpmdiff",                  "osaltpmdiff",                          "ocean_z_salt%4yr_%2mo", "all", "mean", "none", 1
"ocean_model_z",  "boundary_forcing_salt_tendency","boundary_forcing_salt_tendency",      "ocean_z_salt%4yr_%2mo", "all", "mean", "none", 1
"ocean_model_z",  "S_advection_xy",               "S_advection_xy",                       "ocean_z_salt%4yr_%2mo", "all", "mean", "none", 1
"ocean_model_z",  "Sh_tendency_vert_remap",       "Sh_tendency_vert_remap",               "ocean_z_salt%4yr_%2mo", "all", "mean", "none", 1

# SNAP
"ocean_model_z",  "salt",                         "salt"                                  "ocean_z_snap%4yr_%2mo", "all", "none", "none", 2

# ------------------------ #
# Thickness + budget terms #
# ------------------------ #

# MEAN
"ocean_model_z",  "thkcello",                     "thkcello",                             "ocean_z_thck%4yr_%2mo", "all", "mean", "none", 2

"ocean_model_z",  "dhdt",                         "dhdt",                                 "ocean_z_thck%4yr_%2mo", "all", "mean", "none", 1
"ocean_model_z",  "dynamics_h_tendency",          "dynamics_h_tendency",                  "ocean_z_thck%4yr_%2mo", "all", "mean", "none", 1
"ocean_model_z",  "vert_remap_h_tendency",        "vert_remap_h_tendency",                "ocean_z_thck%4yr_%2mo", "all", "mean", "none", 1
"ocean_model_z",  "boundary_forcing_h_tendency",  "boundary_forcing_h_tendency",          "ocean_z_thck%4yr_%2mo", "all", "mean", "none", 1
"ocean_model_z",  "internal_heat_h_tendency",     "internal_heat_h_tendency",             "ocean_z_thck%4yr_%2mo", "all", "mean", "none", 1

# SNAP
"ocean_model_z",  "thkcello",                     "thkcello"                              "ocean_z_snap%4yr_%2mo", "all", "none", "none", 2
"ocean_model_z",  "e",                            "e"                                     "ocean_z_snap%4yr_%2mo", "all", "none", "none", 2

# ----------------------------------------- #
# Other terms necessary for WMT calculation #
# ----------------------------------------- #

# MEAN
"ocean_model_z",  "rhoinsitu",                    "rhoinsitu",                            "ocean_z_thck%4yr_%2mo", "all", "mean", "none", 2
"ocean_model_z",  "rhopot0",                      "rhopot0",                              "ocean_z_thck%4yr_%2mo", "all", "mean", "none", 2
"ocean_model_z",  "rhopot2",                      "rhopot2",                              "ocean_z_thck%4yr_%2mo", "all", "mean", "none", 2
"ocean_model_z",  "drho_dT",                      "drho_dT",                              "ocean_z_thck%4yr_%2mo", "all", "mean", "none", 2
"ocean_model_z",  "drho_dS",                      "drho_dS",                              "ocean_z_thck%4yr_%2mo", "all", "mean", "none", 2
"ocean_model_z",  "umo",                          "umo",                                  "ocean_z_thck%4yr_%2mo", "all", "mean", "none", 2
"ocean_model_z",  "vmo",                          "vmo",                                  "ocean_z_thck%4yr_%2mo", "all", "mean", "none", 2
"ocean_model_z",  "e",                            "e",                                    "ocean_z_thck%4yr_%2mo", "all", "mean", "none", 2

# SNAP
"ocean_model_z",  "rhopot2",                      "rhopot2"                               "ocean_z_snap%4yr_%2mo", "all", "none", "none", 2