The rough plan would be to breakout this section into a new stand-alone function called calc_ts_derivs():
if interfaces is not None:
_ds = xr.Dataset({"thetao": thetao, "so": so})
grid = xgcm.Grid(_ds, coords={"Z": {"center": zcoord}})
thetao = grid.transform(thetao, "Z", interfaces, method="linear")
so = grid.transform(so, "Z", interfaces, method="linear")
zcoord = interfaces.name
pres = (thetao[zcoord] * 1.0e4) + patm
alpha = calc_alpha(thetao, so, pres, eos=eos)
beta = calc_beta(thetao, so, pres, eos=eos)
dtdz = thetao.differentiate(zcoord, edge_order=2)
dsdz = so.differentiate(zcoord, edge_order=2)
SCI will also need mixed layer depth as an input. This would require a completely new function called calc_mld() that can also go in the momlevel.derived module.
Alpha and beta ocean regimes described in Caneill et al. 2022 can be easily implemented by reworking the existing momlevel.derived.calc_n2 function.
The rough plan would be to breakout this section into a new stand-alone function called
calc_ts_derivs():Then,
calc_n2function to call calc_ts_derivsmomlevel.derivedmodule to calculate SCI, which will also callcalc_ts_derivsSCI will also need mixed layer depth as an input. This would require a completely new function called
calc_mld()that can also go in themomlevel.derivedmodule.Things to check:
With: @oceanclimateconnections-temple