An example of a zonal average plot of potential temperature difference between two runs.

demo_zonal_average_T_difference.py

import netCDF4
import numpy
import matplotlib
import matplotlib.pyplot as plt

area = netCDF4.Dataset('/net2/aja/OWG/ocean_static.nc').variables['area_t'][:]
lat = netCDF4.Dataset('/net2/aja/OWG/ocean_static.nc').variables['geolat_c'][:]
T1 = netCDF4.Dataset('/archive/bls/tikal_mom6_2014.01.17/MOM6z_SIS_025_enstro_iTides_1TW/gfdl.ncrc2-intel-prod/pp/ocean_annual/av/annual_20yr/ocean_annual.1900-1919.ann.nc').variables['temp'][0,:]
T2 = netCDF4.Dataset('/archive/bls/tikal_mom6_2014.01.17/MOM6z_SIS_025_enstro_rev5channels/gfdl.ncrc2-intel-prod/pp/ocean_annual/av/annual_20yr/ocean_annual.1900-1919.ann.nc').variables['temp'][0,:]
zInterface = netCDF4.Dataset('/archive/bls/tikal_mom6_2014.01.17/MOM6z_SIS_025_enstro_iTides_1TW/gfdl.ncrc2-intel-prod/pp/ocean_annual/av/annual_20yr/ocean_annual.1900-1919.ann.nc').variables['e'][0,:]

hThickness = zInterface[:-1] - zInterface[1:]
volumes = area * hThickness
delTxave = numpy.sum( volumes * (T1 - T2), axis=-1 ) / numpy.sum( volumes, axis=-1)

y = numpy.amax( lat, axis=-1) # Note: zonal averaging is inappropriate on tri-polar grid
z = numpy.amin( zInterface, axis=-1 )

cmap = plt.get_cmap('seismic')
cmap.set_over([0.3,0.1,0.1])
cmap.set_under([0.1,0.1,0.3])
ci = numpy.linspace(-0.1, 0.1, 21)
norm = matplotlib.colors.BoundaryNorm( ci, cmap.N)
ch = plt.pcolormesh( y, z, delTxave, cmap=cmap, norm=norm );plt.colorbar(ch, extend='both')
plt.clim(-.1,.1)
plt.xlim(-80, 90); plt.ylim(-6500, 0)
plt.xlabel(r'Latitude [$\degree$N]')
plt.ylabel('z* [m]')
plt.title(r'Zonal average $\theta$ difference [$\degree$C]')
plt.show()