<jupyter_start><jupyter_text>Load the dataThe ECCOv4r3 data was converted from its raw MDS (.data / .meta file) format to zarr format, using the [xmitgcm](http://xmitgcm.readthedocs.io) package. [Zarr](http://zarr.readthedocs.io) is a powerful data storage format that can be thought of as an alternative to HDF. In contrast to HDF, zarr works very well with cloud object storage. Zarr is currently useable in python, java, C++, and julia. It is likely that zarr will form the basis of the next major version of the netCDF library.If you're curious, here are some resources to learn more about zarr:  - https://zarr.readthedocs.io/en/stable/tutorial.html  - https://speakerdeck.com/rabernat/pangeo-zarr-cloud-data-storage  - https://mrocklin.github.com/blog/work/2018/02/06/hdf-in-the-cloudThe ECCO zarr data currently lives in [Google Cloud Storage](https://cloud.google.com/storage/) as part of the [Pangeo Data Catalog](http://pangeo.io/catalog.html). This means we can open the whole dataset using one line of code.This takes a bit of time to run because the metadata must be downloaded and parsed. The type of object returned is an [Xarray dataset](http://xarray.pydata.org/en/latest/data-structures.htmldataset).<jupyter_code>cat = intake.open_catalog('catalog.yaml')
list(cat)
ds = cat.ECCOv4r3.to_dask()
ds
ds.nbytes / 1e12
# Load mixed layer depth from my Google cloud bucket 
ds_MXLDEPTH = xr.open_zarr(gcsfs.GCSMap('ecco-data/MXLDEPTH'))<jupyter_output><empty_output><jupyter_text>Note that no data has been actually download yet. Xarray uses the approach of _lazy evaluation_, in which loading of data and execution of computations is delayed as long as possible (i.e. until data is actually needed for a plot). The data are represented symbolically as [dask arrays](http://docs.dask.org/en/latest/array.html). For example:    SALT       (time, k, face, j, i) float32 dask.array    The full shape of the array is (288, 50, 13, 90, 90), quite large. But the _chunksize_ is (1, 50, 13, 90, 90). Here the chunks correspond to the individual granuales of data (objects) in cloud storage. The chunk is the minimum amount of data we can read at one time.<jupyter_code># a trick to make things work a bit faster
coords = ds.coords.to_dataset().reset_coords()
ds = ds.reset_coords(drop=True)<jupyter_output><empty_output><jupyter_text>Visualizing Data A Direct PlotLet's try to visualize something simple: the `Depth` variable. Here is how the data are stored:    Depth      (face, j, i) float32 dask.arrayAlthough depth is a 2D field, there is an extra, dimension (`face`) corresponding to the LLC face number. Let's use xarray's built in plotting functions to plot each face individually.<jupyter_code>coords.Depth.plot(col='face', col_wrap=5)<jupyter_output><empty_output><jupyter_text>This view is not the most useful. It reflects how the data is arranged logically, rather than geographically. A Pretty MapTo make plotting easier, I have created a simple module (`llcmapping.py`, try opening it and looking at the code) that does one thing: regrid the LLC data and plot it on a lat-lon map.<jupyter_code>import importlib
import llcmapping
importlib.reload(llcmapping)
from llcmapping import LLCMapper
#facen=[0,1,2,3,4,5,6,7,8,9,10,11,12]
mapper = LLCMapper(coords)
lons = coords.XC
lats = coords.YC
bnds = [0,359.9,-70,70]
plt.figure(figsize=(12,10))
mapper(coords.Depth, bnds=bnds, cmap='YlGnBu', vmax=6000,vmin=0)
mxldepth = ds_MXLDEPTH.MXLDEPTH
mxldepth.coords['i'] = coords['i']
mxldepth.coords['j'] = coords['j']
mxldepth
# Calculate mean mixed layer depth
mxldepth_clim=mxldepth.mean(dim='time').load()
# Plot mean mixed layer depth
mapper(mxldepth_clim, bnds=bnds, cmap='YlGnBu', vmax=300,vmin=0)<jupyter_output><empty_output><jupyter_text>We can use this with any 2D cell-centered LLC variable. Selecting dataThe entire ECCOv4e3 dataset is contained in a single `Xarray.Dataset` object. How do we find a view specific pieces of data? This is handled by Xarray's [indexing and selecting functions](http://xarray.pydata.org/en/latest/indexing.html). To get the SST from January 2000, we do this:<jupyter_code>sst = ds.THETA.sel(time='2000-01-15', k=0)
sst<jupyter_output><empty_output><jupyter_text>Still no data has been actually downloaded. That doesn't happen until we call `.load()` explicitly or try to make a plot.  Speeding things up with a Dask ClusterHow can we speed things up? In general, the main bottleneck for this type of data analysis is the speed with which we can read the data. With cloud storage, the access is highly parallelizeable.From a Pangeo environment, we can create a [Dask cluster](https://distributed.dask.org/en/latest/) to spread the work out amongst many compute nodes. This works on both HPC and cloud. In the cloud, the compute nodes are provisioned on the fly and can be shut down as soon as we are done with our analysis.The code below will create a cluster with five compute nodes, each with 20GB of RAM.<jupyter_code>#from dask.distributed import Client, progress
#from dask_kubernetes import KubeCluster
#cluster = KubeCluster(n_workers=20)
#cluster
#client = Client(cluster)
#client<jupyter_output><empty_output><jupyter_text>Now we re-run the mean calculation. Note how the dashboard at right helps us visualize what the cluster is doing.<jupyter_code>%time ds.THETA.isel(k=0).mean(dim='time').load()<jupyter_output>CPU times: user 1min, sys: 40.6 s, total: 1min 41s
Wall time: 54 s<jupyter_text>Spatially-Integrated Heat Content AnomalyNow let's do something harder. We will calculate the horizontally integrated heat content anomaly for the full 3D model domain.<jupyter_code># the monthly climatology
theta_clim = ds.THETA.groupby('time.month').mean(dim='time')
# the anomaly
theta_anom = ds.THETA.groupby('time.month') - theta_clim
rho0 = 1029
cp = 3994
ohc = rho0 * cp * (theta_anom *
                   coords.rA *
                   coords.hFacC).sum(dim=['face', 'j', 'i'])
ohc
# actually load the data
ohc.load()
# put the depth coordinate back for plotting purposes
ohc.coords['Z'] = coords.Z
ohc.swap_dims({'k': 'Z'}).transpose().plot(cmap='RdBu_r', vmax=1e20)<jupyter_output><empty_output><jupyter_text>Spatial Derivatives: Heat BudgetAs our final exercise, we will do something much more complicated. We will compute the time-mean convergence of vertically-integrated heat fluxes. This is hard for several reasons.The first reason it is hard is because it involves variables located at different grid points.Following MITgcm conventions, xmitgcm (which produced this dataset) labels the center point with the coordinates `j, i`, the u-velocity point as `j, i_g`, and the v-velocity point as `j_g, i`. The horizontal advective heat flux variables are    ADVx_TH    (time, k, face, j, i_g) float32 dask.array    ADVy_TH    (time, k, face, j_g, i) float32 dask.array    Xarray won't allow us to add or multiply variables that have different dimensions, and xarray by itself doesn't understand how to transform from one grid position to another.**That's why [xgcm](https://xgcm.readthedocs.io/en/latest/) was created.**Xgcm allows us to create a `Grid` object, which understands how to interpolate and take differences in a way that is compatible with finite volume models such at MITgcm. Xgcm also works with many other models, including ROMS, POP, MOM5/6, NEMO, etc.A second reason this is hard is because of the complex topology connecting the different MITgcm faces. Fortunately xgcm also [supports this](https://xgcm.readthedocs.io/en/latest/grid_topology.html).<jupyter_code>import xgcm

# define the connectivity between faces
face_connections = {'face':
                    {0: {'X':  ((12, 'Y', False), (3, 'X', False)),
                         'Y':  (None,             (1, 'Y', False))},
                     1: {'X':  ((11, 'Y', False), (4, 'X', False)),
                         'Y':  ((0, 'Y', False),  (2, 'Y', False))},
                     2: {'X':  ((10, 'Y', False), (5, 'X', False)),
                         'Y':  ((1, 'Y', False),  (6, 'X', False))},
                     3: {'X':  ((0, 'X', False),  (9, 'Y', False)),
                         'Y':  (None,             (4, 'Y', False))},
                     4: {'X':  ((1, 'X', False),  (8, 'Y', False)),
                         'Y':  ((3, 'Y', False),  (5, 'Y', False))},
                     5: {'X':  ((2, 'X', False),  (7, 'Y', False)),
                         'Y':  ((4, 'Y', False),  (6, 'Y', False))},
                     6: {'X':  ((2, 'Y', False),  (7, 'X', False)),
                         'Y':  ((5, 'Y', False),  (10, 'X', False))},
                     7: {'X':  ((6, 'X', False),  (8, 'X', False)),
                         'Y':  ((5, 'X', False),  (10, 'Y', False))},
                     8: {'X':  ((7, 'X', False),  (9, 'X', False)),
                         'Y':  ((4, 'X', False),  (11, 'Y', False))},
                     9: {'X':  ((8, 'X', False),  None),
                         'Y':  ((3, 'X', False),  (12, 'Y', False))},
                     10: {'X': ((6, 'Y', False),  (11, 'X', False)),
                          'Y': ((7, 'Y', False),  (2, 'X', False))},
                     11: {'X': ((10, 'X', False), (12, 'X', False)),
                          'Y': ((8, 'Y', False),  (1, 'X', False))},
                     12: {'X': ((11, 'X', False), None),
                          'Y': ((9, 'Y', False),  (0, 'X', False))}}}

# create the grid object
grid = xgcm.Grid(ds, periodic=False, face_connections=face_connections)
grid<jupyter_output><empty_output><jupyter_text>Now we can use the `grid` object we created to take the divergence of a 2D vector<jupyter_code># vertical integral and time mean of horizontal advective heat flux
advx_th_vint = ds.ADVx_TH.sum(dim='k').mean(dim='time')
advy_th_vint = ds.ADVy_TH.sum(dim='k').mean(dim='time')

# difference in the x and y directions
diff_ADV_th = grid.diff_2d_vector({'X': advx_th_vint, 'Y': advy_th_vint}, boundary='fill')
# convergence
conv_ADV_th = -diff_ADV_th['X'] - diff_ADV_th['Y']
conv_ADV_th
# vertical integral and time mean of horizontal diffusive heat flux
difx_th_vint = ds.DFxE_TH.sum(dim='k').mean(dim='time')
dify_th_vint = ds.DFyE_TH.sum(dim='k').mean(dim='time')

# difference in the x and y directions
diff_DIF_th = grid.diff_2d_vector({'X': difx_th_vint, 'Y': dify_th_vint}, boundary='fill')
# convergence
conv_DIF_th = -diff_DIF_th['X'] - diff_DIF_th['Y']
conv_DIF_th
# convert to Watts / m^2 and load
mean_adv_conv = rho0 * cp * (conv_ADV_th/coords.rA).fillna(0.).load()
mean_dif_conv = rho0 * cp * (conv_DIF_th/coords.rA).fillna(0.).load()
# Standard Deviation of SST
sst_anom = theta_anom.isel(k=0)
sst_std = sst_anom.std(dim='time').load()
# Standard Deviation of Surface heat flux
Qs_clim = ds.TFLUX.groupby('time.month').mean(dim='time')
Qs_anom = ds.TFLUX.groupby('time.month') - Qs_clim
Qs_std = Qs_anom.std(dim='time').load()
# Vertical integral of ocean heat transport over ocean mixed layer
ADVx_TH_clim = ds.ADVx_TH.groupby('time.month').mean(dim='time')
ADVx_TH_anom = ds.ADVx_TH.groupby('time.month') - ADVx_TH_clim

ADVy_TH_clim = ds.ADVy_TH.groupby('time.month').mean(dim='time')
ADVy_TH_anom = ds.ADVy_TH.groupby('time.month') - ADVy_TH_clim

delz = coords.drF*coords.hFacC
v = coords.rA*coords.drF*coords.hFacC

#ADVx_TH_anom = ADVx_TH_anom.where(mxlpoints)
#ADVy_TH_anom = ADVy_TH_anom.where(mxlpoints)

#ADVx_TH_anom_vint = ADVx_TH_anom.sum(dim='k')
#ADVy_TH_anom_vint = ADVy_TH_anom.sum(dim='k')

#delz = coords.drF*coords.hFacC
delz = coords.drF*coords.hFacC

#ADV_TH_diff_anom_vint = grid.diff_2d_vector({'X': ADVx_TH_anom_vint, 'Y': ADVy_TH_anom_vint}, boundary='fill')

ADV_TH_diff_anom = grid.diff_2d_vector({'X': ADVx_TH_anom, 'Y': ADVy_TH_anom}, boundary='fill')

ADV_TH_conv_anom_vint = -((delz/v)*ADV_TH_diff_anom['X']).sum(dim='k') - ((delz/v)*ADV_TH_diff_anom['Y']).sum(dim='k')

ADV_TH_std = ADV_TH_conv_anom_vint.std(dim='time')

ADV_TH_std = rho0 * cp * ADV_TH_std.fillna(0.).load()
#Calculate standard deviation of vertically-integrated ocean diffusive heat convergence
DFxE_TH_clim = ds.DFxE_TH.groupby('time.month').mean(dim='time')
DFxE_TH_anom = ds.DFxE_TH.groupby('time.month') - DFxE_TH_clim

DFyE_TH_clim = ds.DFyE_TH.groupby('time.month').mean(dim='time')
DFyE_TH_anom = ds.DFyE_TH.groupby('time.month') - DFyE_TH_clim

#DFxE_TH_anom_vint = DFxE_TH_anom.sum(dim='k')
#DFyE_TH_anom_vint = DFyE_TH_anom.sum(dim='k')

DF_TH_diff_anom = grid.diff_2d_vector({'X': DFxE_TH_anom, 'Y': DFyE_TH_anom}, boundary='fill')

DF_TH_conv_anom_vint = -((delz/v)*DF_TH_diff_anom['X']).sum(dim='k')-((delz/v)*DF_TH_diff_anom['Y']).sum(dim='k')

DF_TH_std = DF_TH_conv_anom_vint.std(dim='time').load()

DF_TH_std = rho0 * cp * DF_TH_std.fillna(0.).load()

m, ax = mapper(sst_std**2, bnds=bnds, cmap='cubehelix_r', vmax=1.6, vmin=0);
ax.set_title(r'SST Variance (K$^2$)');
m, ax = mapper(mean_adv_conv, bnds=bnds, cmap='RdBu_r', vmax=200, vmin=-200);
ax.set_title(r'Convergence of Advective Flux (W/m$^2$)');
m, ax = mapper(mean_dif_conv,  bnds=bnds, cmap='RdBu_r', vmax=200, vmin=-200)
ax.set_title(r'Convergence of Diffusive Flux (W/m$^2$)');
m, ax = mapper(mean_dif_conv + mean_adv_conv,  bnds=bnds, cmap='RdBu_r', vmax=200, vmin=-200)
ax.set_title(r'Convergence of Net Horizontal Flux (W/m$^2$)');
# Plot standard deviation of vertically-integrated ocean heat convergence
m, ax = mapper(ADV_TH_std,  bnds=bnds, cmap='cubehelix_r', vmax=300, vmin=0);
ax.set_title(r'Standard Deviation of Vertically-Integrated Ocean Heat Convergence (W/m$^2$)');
# Standard deviation of surface heat flux
m, ax = mapper(Qs_std,  bnds=bnds, cmap='cubehelix_r', vmax=60, vmin=0);
ax.set_title(r'Standard Deviation of Surface Heat Flux (W/m$^2$)');
# Standard deviation of vertically-integrated ocean diffusive heat convergence
m, ax = mapper(DF_TH_std,  bnds=bnds, cmap='cubehelix_r', vmax=60, vmin=0);
ax.set_title(r'Standard Deviation of Vertically-Integrated Ocean Diffusive Heat Convergence (W/m$^2$)')
# Test detrending
# import xscale
# from xscale.xscale.signal.fitting import detrend, linreg
import scipy.signal as signal
# Mask 0 entries (continents)
theta_anom = theta_anom.where(theta_anom != 0)
landmask = theta_anom == np.nan
landmask
#Plot global mean time series (not detrended)
areas = coords['rA']
globearea = areas.sum()
weights = areas/globearea
sst_anom = theta_anom.isel(k=0)
sst_globmean = (weights*sst_anom).sum(dim=('face','i','j')).load()
#sst_base = sst_anom.sel(time=slice('1992-01-15','2000-01-15'))
#sst_globmean_base = (weights*sst_base).sum().load()
#sst_globmean_base #Why is this zero? FIX.
plt.figure(figsize=(12,6))
plt.plot(sst_anom.time, sst_globmean)
plt.axhline(0,color='k',linewidth=1)
# Detrend SST
#test = detrend(test, dim='time', type='linear')
#test.dims
sst_anom = sst_anom.fillna(0.)
sst_anom = xr.DataArray(signal.detrend(sst_anom, axis=0), dims=sst_anom.dims, coords=sst_anom.coords)
#Rechunk after detrending
sst_anom.chunk({'time':1})
sst_anom = sst_anom.where(sst_anom != 0)
sst_globmean_detr = sst_anom.mean(dim=('face','i','j')).load()
plt.figure(figsize=(12,6))
plt.plot(sst_anom.time, sst_globmean, label='original')
plt.plot(sst_anom.time, sst_globmean_detr, label='detrended')
plt.legend()
sst_std = sst_anom.std(dim='time')
m, ax = mapper(sst_std,  bnds=bnds, cmap='cubehelix_r', vmax=1.3, vmin=0);
ax.set_title(r'Standard Deviation of Detrended SST (C$^{\circ}$)');
Qs_anom = xr.DataArray(signal.detrend(Qs_anom, axis=0), dims=Qs_anom.dims, coords=Qs_anom.coords)
#Rechunk after detrending
Qs_anom.chunk({'time':1})
Qs_std_detr = Qs_anom.std(dim='time').load()
# Standard deviation of surface heat flux
m, ax = mapper(Qs_std_detr,  bnds=bnds, cmap='cubehelix_r', vmax=60, vmin=0);
ax.set_title(r'Standard Deviation of Detrended Surface Heat Flux (W/m$^2$)');
# Test subsetting ocean mixed layer points
mxlpoints = np.abs(coords['Z']) <= mxldepth_clim
#delz = coords.drF*coords.hFacC
delz = coords.drF*coords.hFacC
delz=delz.where(mxlpoints)
delz_sum = delz.sum(dim='k')
mxl_theta = theta_anom.where(mxlpoints)
mxl_theta
weights = delz/delz_sum
#mxl_theta_mean = mxl_theta.mean(dim='k', skipna=True)
mxl_theta_mean = (mxl_theta*weights).sum(dim='k')
mxl_theta_mean
mxl_theta_std = mxl_theta_mean.std(dim='time').load()
m, ax = mapper(mxl_theta_std,  bnds=bnds, cmap='cubehelix_r', vmax=1.3, vmin=0);
ax.set_title(r'Variance of Mean Mixed Layer Temperature (C$^{\circ}$)');
m, ax = mapper(sst_std,  bnds=bnds, cmap='cubehelix_r', vmax=1.3, vmin=0);
ax.set_title(r' SST (C$^{\circ}$)');
# Mixed layer vertical integral of time-averaged vertical heat transport
v = coords.rA*coords.drF*coords.hFacC
delz = coords.drF*coords.hFacC
ADVr_TH_clim = ds.ADVr_TH.mean(dim='time')
ADVr_TH_diff_clim = grid.diff(ADVr_TH_clim, 'Z', boundary='extend')
ADVr_TH_clim_vint = ((delz/v)*ADVr_TH_diff_clim).where(mxlpoints).sum(dim='k')

ADVr_TH_mean = rho0 * cp * ADVr_TH_clim_vint.fillna(0.).load()
# Mixed layer vertical integral of time-averaged horizontal heat transport
ADVx_TH_clim = ds.ADVx_TH.mean(dim='time')
ADVy_TH_clim = ds.ADVy_TH.mean(dim='time')
ADV_TH_diff_clim = grid.diff_2d_vector({'X': ADVx_TH_clim, 'Y': ADVy_TH_clim}, boundary='fill')
ADV_TH_conv_clim_vint = -((delz/v)*ADV_TH_diff_clim['X']).where(mxlpoints).sum(dim='k') - ((delz/v)*ADV_TH_diff_clim['Y']).where(mxlpoints).sum(dim='k')
ADV_TH_mean = rho0 * cp * ADV_TH_conv_clim_vint.fillna(0.).load()
# Same as previous except for horizontal diffusion 
DFxE_TH_clim = ds.DFxE_TH.mean(dim='time')
DFyE_TH_clim = ds.DFyE_TH.mean(dim='time')
DF_TH_diff_clim = grid.diff_2d_vector({'X': DFxE_TH_clim, 'Y': DFyE_TH_clim}, boundary='fill')
DF_TH_conv_clim_vint = -((delz/v)*DF_TH_diff_clim['X']).where(mxlpoints).sum(dim='k') - ((delz/v)*DF_TH_diff_clim['Y']).where(mxlpoints).sum(dim='k')
DF_TH_mean = rho0 * cp * DF_TH_conv_clim_vint.fillna(0.).load()


# Same as previous except for vertical diffusion (explicit and implicit)
DFrE_TH_clim = ds.DFrE_TH.mean(dim='time')
DFrI_TH_clim = ds.DFrI_TH.mean(dim='time')
DFrE_TH_diff_clim = grid.diff(DFrE_TH_clim, 'Z', boundary='extend')
DFrI_TH_diff_clim = grid.diff(DFrI_TH_clim, 'Z', boundary='extend')

DFrE_TH_clim_vint = ((delz/v)*DFrE_TH_diff_clim).where(mxlpoints).sum(dim='k')
DFrI_TH_clim_vint = ((delz/v)*DFrI_TH_diff_clim).where(mxlpoints).sum(dim='k')

DFr_TH_clim_vint = DFrE_TH_clim_vint + DFrI_TH_clim_vint

DFr_TH_mean = rho0 * cp * DFr_TH_clim_vint.fillna(0.).load()
# Mixed layer vertical integral of ocean heat transport 
ADVx_TH_clim = ds.ADVx_TH.groupby('time.month').mean(dim='time')
ADVx_TH_anom = ds.ADVx_TH.groupby('time.month') - ADVx_TH_clim

ADVy_TH_clim = ds.ADVy_TH.groupby('time.month').mean(dim='time')
ADVy_TH_anom = ds.ADVy_TH.groupby('time.month') - ADVy_TH_clim

ADV_TH_diff_anom = grid.diff_2d_vector({'X': ADVx_TH_anom, 'Y': ADVy_TH_anom}, boundary='fill')

ADV_TH_conv_anom_vint = -((delz/v)*ADV_TH_diff_anom['X']).where(mxlpoints).sum(dim='k') - ((delz/v)*ADV_TH_diff_anom['Y']).where(mxlpoints).sum(dim='k')

ADV_TH_std = ADV_TH_conv_anom_vint.std(dim='time')

ADV_TH_std = rho0 * cp * ADV_TH_std.fillna(0.).load()
# Mixed layer vertical integral of vertical ocean heat convergence 
ADVr_TH_clim = ds.ADVr_TH.groupby('time.month').mean(dim='time')
ADVr_TH_anom = ds.ADVr_TH.groupby('time.month') - ADVr_TH_clim

ADVr_TH_diff_anom = grid.diff(ADVr_TH_anom, 'Z', boundary='extend')

ADVr_TH_conv_anom_vint = ((delz/v)*ADVr_TH_diff_anom).where(mxlpoints).sum(dim='k')

ADVr_TH_std = ADVr_TH_conv_anom_vint.std(dim='time')

ADVr_TH_std = rho0 * cp * ADVr_TH_std.fillna(0.).load()
# Same as previous except for diffusive terms
DFxE_TH_clim = ds.DFxE_TH.groupby('time.month').mean(dim='time')
DFxE_TH_anom = ds.DFxE_TH.groupby('time.month') - DFxE_TH_clim

DFyE_TH_clim = ds.DFyE_TH.groupby('time.month').mean(dim='time')
DFyE_TH_anom = ds.DFyE_TH.groupby('time.month') - DFyE_TH_clim

DF_TH_diff_anom = grid.diff_2d_vector({'X': DFxE_TH_anom, 'Y': DFyE_TH_anom}, boundary='fill')

DF_TH_conv_anom_vint = -((delz/v)*DF_TH_diff_anom['X']).where(mxlpoints).sum(dim='k') - ((delz/v)*DF_TH_diff_anom['Y']).where(mxlpoints).sum(dim='k')

DF_TH_std = DF_TH_conv_anom_vint.std(dim='time')

DF_TH_std = rho0 * cp * DF_TH_std.fillna(0.).load()
DFrE_TH_clim = ds.DFrE_TH.groupby('time.month').mean(dim='time')
DFrE_TH_anom = ds.DFrE_TH.groupby('time.month') - DFrE_TH_clim

DFrI_TH_clim = ds.DFrI_TH.groupby('time.month').mean(dim='time')
DFrI_TH_anom = ds.DFrI_TH.groupby('time.month') - DFrI_TH_clim

DFrE_TH_diff_anom = grid.diff(DFrE_TH_anom, 'Z', boundary='extend')
DFrI_TH_diff_anom = grid.diff(DFrI_TH_anom, 'Z', boundary='extend')

DFrE_TH_conv_anom_vint = ((delz/v)*DFrE_TH_diff_anom).where(mxlpoints).sum(dim='k')
DFrI_TH_conv_anom_vint = ((delz/v)*DFrI_TH_diff_anom).where(mxlpoints).sum(dim='k')

DFr_TH_conv_anom_vint = DFrE_TH_conv_anom_vint + DFrI_TH_conv_anom_vint

DFr_TH_std = DFr_TH_conv_anom_vint.std(dim='time')

DFr_TH_std = rho0 * cp * DFr_TH_std.fillna(0.).load()

# Time-averaged vertical heat convergence in ocean mixed layer
m, ax = mapper(ADVr_TH_mean,  bnds=bnds, cmap='RdBu_r', vmax=500, vmin=-500);
ax.set_title(r'Time-averaged Vertical Heat Convergence in Mixed Layer (W/m$^2$)');
# Time-averaged horizontal heat convergence in ocean mixed layer
m, ax = mapper(ADV_TH_mean,  bnds=bnds, cmap='RdBu_r', vmax=500, vmin=-500);
ax.set_title(r'Time-averaged Horizontal Heat Convergence in Mixed Layer (W/m$^2$)');
# Time-averaged horizontal diffusive heat convergence in ocean mixed layer
m, ax = mapper(DF_TH_mean,  bnds=bnds, cmap='RdBu_r', vmax=500, vmin=-500);
ax.set_title(r'Time-averaged Horizontal Diffusive Heat Convergence in Mixed Layer (W/m$^2$)');
# Time-averaged vertical diffusive heat convergence in ocean mixed layer
m, ax = mapper(DFr_TH_mean,  bnds=bnds, cmap='RdBu_r', vmax=500, vmin=-500);
ax.set_title(r'Time-averaged Vertical Diffusive Heat Convergence in Mixed Layer (W/m$^2$)');
# Plot standard deviation of ocean mixed layer horizontal heat convergence
m, ax = mapper(ADV_TH_std,  bnds=bnds, cmap='cubehelix_r', vmax=600, vmin=0);
ax.set_title(r'Standard Deviation of Mixed Layer Horizontal Heat Convergence (W/m$^2$)');
# Plot standard deviation of ocean mixed layer vertical heat convergence
m, ax = mapper(ADVr_TH_std,  bnds=bnds, cmap='cubehelix_r', vmax=600, vmin=0);
ax.set_title(r'Standard Deviation of Ocean Mixed Layer Vertical Heat Convergence (W/m$^2$)');
# Plot standard deviation of horizontal diffusive heat convergence
m, ax = mapper(DF_TH_std,  bnds=bnds, cmap='cubehelix_r', vmax=600, vmin=0);
ax.set_title(r'Standard Deviation of Ocean Mixed Layer Horizontal Diffusive Heat Convergence (W/m$^2$)');
# Plot standard deviation of horizontal diffusive heat convergence
m, ax = mapper(DFr_TH_std,  bnds=bnds, cmap='cubehelix_r', vmax=600, vmin=0);
ax.set_title(r'Standard Deviation of Ocean Mixed Layer Vertical Diffusive Heat Convergence (W/m$^2$)');
# Plot standard deviation of ocean mixed layer horizontal heat convergence
m, ax = mapper(Qs_std,  bnds=bnds, cmap='cubehelix_r', vmax=600, vmin=0);
ax.set_title(r'Standard Deviation of Surface Heat Flux (W/m$^2$)');
# Test Butterworth filtering
# Check energy budget for a grid point (follow psuedocode)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14652978475607126, 'fraction_numerical': 0.016305051115041193, 'mean_word_length': 5.08983451536643, 'pattern_counts': {'":': 0, '<': 41, '<?xml version=': 0, '>': 40, 'https://': 7, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '200587', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 9447, 'n_tokens_neox': 9041, 'n_words': 2097, 'path': 'ecco_meanfields.ipynb', 'repo_name': 'cpatrizio88/pangeo_binder_example', 'size': '23183'}