In [1]:
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import matplotlib.colors as colors
from matplotlib import gridspec
from brokenaxes import brokenaxes
import pandas as pd
import copy
In [2]:
df = pd.read_csv('drifts.csv')
df.set_index(df['model'] + ' (' + df['run'] + ')', drop=True, inplace=True)
df
Out[2]:
In [3]:
def plot_abline(ax, slope, intercept, static_bounds=True):
"""Plot a line from slope and intercept"""
xlim = ax.get_xlim()
ylim = ax.get_ylim()
if type(xlim[0]) in (list, tuple):
for lims in xlim:
x_vals = np.array(lims)
y_vals = intercept + slope * x_vals
ax.plot(x_vals, y_vals, linestyle='--', c='0.5')
else:
x_vals = np.array(xlim)
y_vals = intercept + slope * x_vals
ax.plot(x_vals, y_vals, linestyle='--', c='0.5')
if static_bounds:
ax.set_xlim(xlim)
ax.set_ylim(ylim)
def plot_shading(ax):
"""Plot shading to indicate dominant source of drift."""
xlim = ax.get_xlim()
ylim = ax.get_ylim()
x_vals = np.array(xlim)
y_vals = x_vals * 2
ax.fill_between(x_vals, 0, y_vals, alpha=0.3, color='0.5')
ax.set_xlim(xlim)
ax.set_ylim(ylim)
def plot_eei_shading(ax):
"""Plot shading to indicate netTOA / OHC valid range."""
xlim = ax.get_xlim()
ylim = ax.get_ylim()
x_vals = np.array(xlim)
y_vals = x_vals * 0.8
ax.fill_between(x_vals, x_vals, y_vals, alpha=0.3, color='0.5')
ax.set_xlim(xlim)
ax.set_ylim(ylim)
def format_axis_label(orig_label, units, scale_factor):
"""Put LaTeX math into axis labels"""
label = orig_label.split('(')[0] + '(' + units + ')'
label = label.replace('(', '($').replace(')', '$)')
label = label.replace('s-1', '\; s^{-1}')
label = label.replace('m-2', '\; m^{-2}')
label = label.replace('yr-1', '\; yr^{-1}')
if scale_factor:
scale_factor = int(scale_factor) * -1
label = label.replace('($', '($10^{%s} \;' %(str(scale_factor)))
return label
In [11]:
zoom_limits = {'thermal energy conservation': [-0.15, 0.15],
'mass conservation': [-2e7, 2e7],
'salt conservation': [-7e-13, 7e-13],
'barystatic consistency': [-1e-12, 1e-12],
'planetary energy imbalance': [-0.15, 0.15]}
ocean_model_colors = {'MOM': 'red',
'GOLD': 'gold',
'NEMO': 'blue',
'OPA': 'green',
'COCO': 'chocolate',
'MPI-OM': 'purple',
'MICOM-HAMOCC': 'teal',
'POP': 'lime'}
markers = ['o', '<', '^', '>', 'v', 's', 'p', 'D',
'o', '<', '^', '>', 'v', 's', 'p', 'D',
'o', '<', '^', '>', 'v', 's', 'p', 'D']
def plot_aesthetics(ax, yvar, xvar, units, scinotation, shading, scale_factor,
xpad=None, ypad=None, non_square=True):
"""Set the plot aesthetics"""
if shading:
plot_shading(ax)
if 'netTOA' in xvar:
plot_eei_shading(ax)
else:
plot_abline(ax, 1, 0, static_bounds=non_square)
ax.axhline(y=0, color='0.5', linewidth=1.0)
ax.axvline(x=0, color='0.5', linewidth=1.0)
#ax.yaxis.major.formatter._useMathText = True
#ax.xaxis.major.formatter._useMathText = True
ylabel = format_axis_label(yvar, units, scale_factor)
if ypad:
ax.set_ylabel(ylabel, labelpad=ypad)
else:
ax.set_ylabel(ylabel)
xlabel = format_axis_label(xvar, units, scale_factor)
if xpad:
ax.set_xlabel(xlabel, labelpad=xpad)
else:
ax.set_xlabel(xlabel)
ax.set_xlabel(xlabel, labelpad=xpad)
#plt.sca(ax)
if scinotation:
plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0), useMathText=True)
plt.ticklabel_format(style='sci', axis='x', scilimits=(0,0), useMathText=True)
# Shrink current axis by 20%
#box = ax.get_position()
#ax.set_position([box.x0, box.y0, box.width * 0.8, box.height])
def get_units(column_header):
"""Get the units from the column header."""
units = column_header.split('(')[-1].split(')')[0]
return units
def convert_units(value, start_units, end_units, ocean_area=None):
"""Convert units."""
sec_in_year = 365.25 * 24 * 60 * 60
if start_units == end_units:
new_value = value
else:
assert start_units in ['J yr-1', 'm yr-1', 'kg yr-1', 'g/kg yr-1', 'm yr-1']
assert end_units in ['PW', 'W m-2', 'mm yr-1', 'kg s-1', 'g/kg s-1', 'm s-1']
if start_units == 'J yr-1':
new_value = value / sec_in_year
if end_units == 'W m-2':
earth_surface_area = 5.1e14
new_value = new_value / earth_surface_area
elif end_units == 'PW':
new_value = new_value / 1e15
elif (start_units == 'm yr-1') and (end_units == 'mm yr-1'):
new_value = value * 1000
elif (start_units == 'kg yr-1') and (end_units == 'mm yr-1'):
assert ocean_area
new_value = value / ocean_area
elif (start_units == 'kg yr-1') and (end_units == 'kg s-1'):
new_value = value / sec_in_year
elif (start_units == 'g/kg yr-1') and (end_units == 'g/kg s-1'):
new_value = value / sec_in_year
elif (start_units == 'm yr-1') and (end_units == 'm s-1'):
new_value = value / sec_in_year
return new_value
def plot_comparison(df, title, xvar, yvar, plot_units, scale_factor=0,
scinotation=False, shading=False, zoom=None, outfile=None):
"""Plot comparison for given x and y variables.
Data are multiplied by 10^scale_factor.
"""
if zoom:
fig = plt.figure(figsize=[14, 5])
ax1 = fig.add_subplot(1, 2, 1)
ax2 = fig.add_subplot(1, 2, 2)
else:
fig = plt.figure(figsize=[7, 5])
ax1 = fig.add_subplot(1, 1, 1)
#colormap = cm.gist_rainbow
#colorlist = [colors.rgb2hex(colormap(i)) for i in np.linspace(0, 0.9, len(df['model']))]
x_input_units = get_units(xvar)
y_input_units = get_units(yvar)
for dotnum in range(len(df['model'])):
area = df['ocean area (m2)'][dotnum]
x = convert_units(df[xvar][dotnum], x_input_units, plot_units, ocean_area=area) * 10**scale_factor
y = convert_units(df[yvar][dotnum], y_input_units, plot_units, ocean_area=area) * 10**scale_factor
marker = markers[dotnum]
label = df['model'][dotnum] + ' (' + df['run'][dotnum] + ')'
ocean_model = df['ocean model'][dotnum]
color = ocean_model_colors[ocean_model]
if df['project'][dotnum] == 'cmip6':
facecolors = color
edgecolors ='none'
else:
facecolors = 'none'
edgecolors = color
ax1.scatter(x, y, label=label, s=130, linewidth=1.2, marker=marker,
facecolors=facecolors, edgecolors=edgecolors)
if zoom:
ax2.scatter(x, y, label=label, s=130, linewidth=1.2, marker=marker,
facecolors=facecolors, edgecolors=edgecolors)
plot_aesthetics(ax1, yvar, xvar, plot_units, scinotation, shading, scale_factor)
if zoom:
ax2.set_xlim(zoom)
ax2.set_ylim(zoom)
plot_aesthetics(ax2, yvar, xvar, plot_units, scinotation, shading, scale_factor)
ax1.set_title('all models')
ax2.set_title('zoomed in')
plt.suptitle(title)
ax2.legend(loc='center left', bbox_to_anchor=(1, 0.5))
else:
plt.title(title)
ax1.legend(loc='center left', bbox_to_anchor=(1, 0.5))
if outfile:
plt.savefig(outfile, bbox_inches='tight', dpi=400)
plt.show()
else:
plt.show()
def plot_regression_coefficients(df, comparison_list, title, outfile=None, ylim=None):
"""Plot the regression coefficient for each model"""
labels = {"netTOA vs hfds regression": "cumulative netTOA ($Q_r$) vs. cumulative surface heat flux ($Q_h$)",
"hfds vs thermal OHC regression": "cumulative surface heat flux ($Q_h$) vs. thermal OHC ($H_T$)",
"netTOA vs thermal OHC regression": "cumulative netTOA ($Q_r$) vs. thermal OHC ($H_T$)",
"wfo vs masso regression": "cumulative freshwater flux ($Q_m$) vs. ocean mass ($M$)",
"wfo vs soga regression": "cumulative freshwater flux ($Q_m$) vs. ocean salinity ($S$)",
"masso vs soga regression": "ocean mass ($M$) vs. ocean salinity ($S$)"}
nvars = len(comparison_list)
assert nvars <= 3
x = np.arange(df.shape[0])
if nvars in [1, 3]:
x0 = x - 0.2
x1 = x + 0.2
xlist = [x, x0, x1]
else:
x0 = x - 0.1
x1 = x + 0.1
xlist = [x0, x1]
plt.figure(figsize=[14, 5])
plt.axhline(y=1.0, color='0.5', linewidth=0.5)
plt.axvline(x=14.5, color='0.5', linewidth=2.0)
colors = ['cadetblue', 'seagreen', 'mediumpurple']
linestyles = ['--', '-.', ':']
for pnum, var in enumerate(comparison_list):
y = df[var].to_numpy()
color = colors[pnum]
linestyle = linestyles[pnum]
markerline, stemlines, baseline = plt.stem(xlist[pnum], y, color, basefmt=" ", bottom=1.0, label=labels[var])
plt.setp(stemlines, color=color, linestyle=linestyle)
plt.setp(markerline, color=color)
if ylim:
plt.ylim(ylim)
plt.axvline(x=x[0]-0.5, color='0.5', linewidth=0.5)
for val in x:
plt.axvline(x=val+0.5, color='0.5', linewidth=0.5)
#plt.grid(True, axis='x')
plt.legend()
plt.xticks(x, df.index.to_list(), rotation=90)
plt.ylabel('regression coefficient')
plt.title(title)
if outfile:
plt.savefig(outfile, bbox_inches='tight', dpi=400)
plt.show()
else:
plt.show()
In [5]:
def plot_broken_comparison(df, title, xvar, yvar, plot_units, scale_factor=0,
scinotation=False, shading=False, outfile=None,
xlims=None, ylims=None, xpad=None, ypad=None,
hspace=0.04, wspace=0.04):
"""Plot comparison for given x and y variables.
Data are multiplied by 10^scale_factor.
"""
fig = plt.figure(figsize=[10, 8])
if xlims and ylims:
bax = brokenaxes(xlims=xlims, ylims=ylims, hspace=hspace, wspace=wspace)
else:
bax = fig.add_subplot(1, 1, 1)
x_input_units = get_units(xvar)
y_input_units = get_units(yvar)
for dotnum in range(len(df['model'])):
area = df['ocean area (m2)'][dotnum]
x = convert_units(df[xvar][dotnum], x_input_units, plot_units, ocean_area=area) * 10**scale_factor
y = convert_units(df[yvar][dotnum], y_input_units, plot_units, ocean_area=area) * 10**scale_factor
marker = markers[dotnum]
label = df['model'][dotnum] + ' (' + df['run'][dotnum] + ')'
ocean_model = df['ocean model'][dotnum]
color = ocean_model_colors[ocean_model]
if df['project'][dotnum] == 'cmip6':
facecolors = color
edgecolors ='none'
else:
facecolors = 'none'
edgecolors = color
bax.scatter(x, y, label=label, s=130, linewidth=1.2, marker=marker,
facecolors=facecolors, edgecolors=edgecolors)
print(label, ((x-y) / 1.8) * 100, '%')
if xlims:
non_square = False
else:
non_square = True
bax.spines["top"].set_visible(False)
bax.spines["right"].set_visible(False)
plot_aesthetics(bax, yvar, xvar, plot_units, scinotation, shading, scale_factor,
xpad=xpad, ypad=ypad, non_square=non_square)
plt.title(title)
bax.legend(loc='center left', bbox_to_anchor=(1, 0.5))
if outfile:
plt.savefig(outfile, bbox_inches='tight', dpi=400)
plt.show()
else:
plt.show()
In [158]:
plot_broken_comparison(df, 'thermal energy conservation', 'hfds (J yr-1)', 'thermal OHC (J yr-1)', 'W m-2')
In [159]:
plot_broken_comparison(df, 'thermal energy conservation', 'hfds (J yr-1)', 'thermal OHC (J yr-1)', 'W m-2',
xlims=[(-41.05, -40.82), (-0.35, 0.3)], ylims=[(-0.3, 0.25), (0.55, 0.65)],
wspace=0.08, hspace=0.08, xpad=20,
outfile='thermal_conservation.png')
In [141]:
plot_comparison(df, 'thermal energy conservation', 'hfds (J yr-1)', 'thermal OHC (J yr-1)', 'PW',
shading=False, zoom=[-0.2, 0.2])
Below the 1:1 line indicates that the ocean model is losing heat.
Within the shaded area means the drift in surface forcing is greater than the internally generated ocean model drift.
In [14]:
plot_comparison(df, 'thermal energy conservation', 'hfds (J yr-1)', 'thermal OHC (J yr-1)', 'W m-2',
shading=False, zoom=[-0.35, 0.35])
The planetary energy imbalance associated with climate change is 0.5 - 1.0 $W m^{-2}$.
In [145]:
plot_regression_coefficients(df, ['hfds vs thermal OHC regression'],
'thermal energy conservation after drift removal',
outfile='thermal_regression.png')
In [10]:
plot_broken_comparison(df, 'mass conservation', 'wfo (kg yr-1)', 'masso (kg yr-1)', 'mm yr-1',
outfile='mass_conservation.png')
In [147]:
plot_comparison(df, 'mass conservation', 'wfo (kg yr-1)', 'masso (kg yr-1)', 'kg s-1',
scinotation=True, shading=False, zoom=[-1e7, 1e7])
Below the 1:1 line indicates that the ocean model is losing mass.
Within the shaded area means the drift in surface forcing is greater than the internally generated ocean model drift.
In [17]:
plot_comparison(df, 'mass conservation', 'wfo (kg yr-1)', 'masso (kg yr-1)', 'mm yr-1',
shading=False, zoom=[-0.7, 0.7])
The current rate of global mean sea level rise is 3.4 mm/year.
An alternative unit for comparison (which has the advantage of avoiding using the global ocean area, which is different for different models) might be Gt/year (to compare with Antarctic and/or Greenland melt rates).
In [55]:
plot_regression_coefficients(df, ['wfo vs masso regression'],
'mass conservation after drift removal')
In [132]:
plot_broken_comparison(df, 'salt conservation', 'masso (g/kg yr-1)', 'soga (g/kg yr-1)', 'g/kg s-1', scale_factor=13)
In [164]:
plot_broken_comparison(df, 'salt conservation', 'masso (g/kg yr-1)', 'soga (g/kg yr-1)', 'g/kg s-1', scale_factor=13,
xlims=[(-2, 5)], ylims=[(-19, -17.5), (-2.2, 5)], hspace=0.1, outfile='salt_conservation.png')
In [19]:
plot_comparison(df, 'salt conservation', 'masso (g/kg yr-1)', 'soga (g/kg yr-1)', 'g/kg s-1',
scale_factor=13, zoom=[-1.2, 1.2])
In [8]:
def plot_broken_regression_coefficients(df, comparison_list, title, outfile=None, ylim=None):
"""Plot the regression coefficient for each model"""
labels = {"netTOA vs hfds regression": "cumulative netTOA ($Q_r$) vs. cumulative surface heat flux ($Q_h$)",
"hfds vs thermal OHC regression": "cumulative surface heat flux ($Q_h$) vs. thermal OHC ($H_T$)",
"netTOA vs thermal OHC regression": "cumulative netTOA ($Q_r$) vs. thermal OHC ($H_T$)",
"wfo vs masso regression": "cumulative freshwater flux ($Q_m$) vs. ocean mass ($M$)",
"wfo vs soga regression": "cumulative freshwater flux ($Q_m$) vs. ocean salinity ($S$)",
"masso vs soga regression": "ocean mass ($M$) vs. ocean salinity ($S$)"}
nvars = len(comparison_list)
assert nvars <= 3
x = np.arange(df.shape[0])
if nvars in [1, 3]:
x0 = x - 0.2
x1 = x + 0.2
xlist = [x, x0, x1]
else:
x0 = x - 0.1
x1 = x + 0.1
xlist = [x0, x1]
fig, (ax1, ax2, ax3) = plt.subplots(3, 1, sharex=True, figsize=(18,8), gridspec_kw={'height_ratios': [1, 5, 1]})
ax1.spines['bottom'].set_visible(False)
ax2.spines['bottom'].set_visible(False)
ax1.tick_params(axis='x', which='both', bottom=False)
ax2.tick_params(axis='x', which='both', bottom=False)
ax2.spines['top'].set_visible(False)
ax3.spines['top'].set_visible(False)
ax1.set_ylim(1.5, 3.0)
ax2.set_ylim(-0.1, 1.4)
ax3.set_ylim(-5.5, -1)
ax2.axhline(y=1.0, color='0.5', linewidth=0.5)
ax1.axvline(x=14.5, color='0.5', linewidth=2.0)
ax2.axvline(x=14.5, color='0.5', linewidth=2.0)
ax3.axvline(x=14.5, color='0.5', linewidth=2.0)
colors = ['cadetblue', 'seagreen', 'mediumpurple']
linestyles = ['--', '-.', ':']
for pnum, var in enumerate(comparison_list):
y = df[var].to_numpy()
color = colors[pnum]
linestyle = linestyles[pnum]
markerline, stemlines, baseline = ax2.stem(xlist[pnum], y, color, basefmt=" ", bottom=1.0, label=labels[var])
plt.setp(stemlines, color=color, linestyle=linestyle)
plt.setp(markerline, color=color)
markerline, stemlines, baseline = ax1.stem(xlist[pnum], y, color, basefmt=" ", bottom=1.0, label=labels[var])
plt.setp(stemlines, color=color, linestyle=linestyle)
plt.setp(markerline, color=color)
markerline, stemlines, baseline = ax3.stem(xlist[pnum], y, color, basefmt=" ", bottom=1.0, label=labels[var])
plt.setp(stemlines, color=color, linestyle=linestyle)
plt.setp(markerline, color=color)
ax1.axvline(x=x[0]-0.5, color='0.5', linewidth=0.5)
ax2.axvline(x=x[0]-0.5, color='0.5', linewidth=0.5)
ax3.axvline(x=x[0]-0.5, color='0.5', linewidth=0.5)
for val in x:
ax1.axvline(x=val+0.5, color='0.5', linewidth=0.5)
ax2.axvline(x=val+0.5, color='0.5', linewidth=0.5)
ax3.axvline(x=val+0.5, color='0.5', linewidth=0.5)
#plt.grid(True, axis='x')
ax1.legend()
plt.xticks(x, df.index.to_list(), rotation=90)
ax2.set_ylabel('regression coefficient')
ax1.set_title(title)
plt.subplots_adjust(hspace=0.1)
if outfile:
plt.savefig(outfile, bbox_inches='tight', dpi=400)
plt.show()
else:
plt.show()
In [9]:
plot_broken_regression_coefficients(df, ['wfo vs masso regression', 'wfo vs soga regression', 'masso vs soga regression'],
'mass and salt conservation after drift removal',
outfile='mass_regression.png')
In [6]:
plot_regression_coefficients(df, ['wfo vs masso regression', 'wfo vs soga regression', 'masso vs soga regression'],
'mass and salt conservation after drift removal',
outfile='mass_regression.png')
#ylim=[-0.2, 1.2]
In [21]:
plot_comparison(df, 'planetary energy imbalance', 'netTOA (J yr-1)', 'thermal OHC (J yr-1)', 'PW')
Grey shading shows the region where $\delta H / \delta t$ = 80-100% of netTOA. Below that regions indicated that the model is losing energy between the top of the atmosphere and ocean (most models lose heat).
In [168]:
plot_broken_comparison(df, 'planetary energy imbalance', 'netTOA (J yr-1)', 'thermal OHC (J yr-1)', 'W m-2',
outfile='eei_conservation.png')
In [12]:
plot_regression_coefficients(df, ['netTOA vs hfds regression', 'hfds vs thermal OHC regression', 'netTOA vs thermal OHC regression'],
'thermal energy conservation after drift removal', outfile='eei_regression.png')
In [16]:
plot_comparison(df, 'barystatic consistency', 'masso (m yr-1)', 'zosbary (m yr-1)', 'm s-1',
scale_factor=11, zoom=[-0.12, 0.12])
In [13]:
df_ohc = df[['OHC (J yr-1)', 'thermal OHC (J yr-1)', 'barystatic OHC (J yr-1)']]
df_ohc.set_index(df['model'] + ' (' + df['run'] + ')', drop=True, inplace=True)
In [14]:
sec_in_year = 365.25 * 24 * 60 * 60
earth_surface_area = 5.1e14
df_ohc = (df_ohc / sec_in_year) / earth_surface_area
df_ohc = df_ohc.rename(columns={"OHC (J yr-1)": "change in OHC ($dH/dt$)",
"thermal OHC (J yr-1)": "change in thermal OHC ($dH_T/dt$)",
"barystatic OHC (J yr-1)": "change in barystatic OHC ($dH_m/dt$)"})
In [15]:
df_ohc.plot.bar(figsize=(18,6), color=['black', 'red', 'blue'], width=0.9)
plt.axvline(x=14.5, color='0.5', linewidth=2.0)
plt.axhline(y=0.5, color='0.5', linewidth=0.5, linestyle='--')
#plt.title('drift in OHC and its components')
plt.ylabel('$W \; m^{-2}$')
plt.savefig('ohc_drift_breakdown.png', bbox_inches='tight', dpi=400)
#plt.show()
In [181]:
(df_ohc['thermal OHC'] / df_ohc['OHC'] ) * 100
Out[181]:
In [16]:
df['atmos energy leakage (J yr-1)'] = df['netTOA (J yr-1)'] - df['hfds (J yr-1)']
df['ocean energy leakage (J yr-1)'] = df['hfds (J yr-1)'] - df['thermal OHC (J yr-1)']
df['total energy leakage (J yr-1)'] = df['netTOA (J yr-1)'] - df['thermal OHC (J yr-1)']
In [17]:
df_leakage = df[['total energy leakage (J yr-1)', 'atmos energy leakage (J yr-1)', 'ocean energy leakage (J yr-1)']]
df_leakage.set_index(df['model'] + ' (' + df['run'] + ')', drop=True, inplace=True)
In [18]:
sec_in_year = 365.25 * 24 * 60 * 60
earth_surface_area = 5.1e14
df_leakage = (df_leakage / sec_in_year) / earth_surface_area
df_leakage = df_leakage.rename(columns={"total energy leakage (J yr-1)": "total leakage ($dQ_r/dt - dH_T/dt$)",
"ocean energy leakage (J yr-1)": "ocean leakage ($dQ_h/dt - dH_T/dt$)",
"atmos energy leakage (J yr-1)": "non-ocean leakage ($dQ_r/dt - dQ_h/dt$)",
})
In [19]:
fig, (ax1, ax2, ax3) = plt.subplots(3, 1, sharex=True, figsize=(18,8), gridspec_kw={'height_ratios': [1, 5, 1]})
ax1.spines['bottom'].set_visible(False)
ax2.spines['bottom'].set_visible(False)
ax1.tick_params(axis='x', which='both', bottom=False)
ax2.tick_params(axis='x', which='both', bottom=False)
ax2.spines['top'].set_visible(False)
ax3.spines['top'].set_visible(False)
ax1.set_ylim(30, 45)
ax2.set_ylim(-4, 3)
ax3.set_ylim(-45, -30)
df_leakage.plot(ax=ax3, kind='bar', color=['gold', 'green', 'blue'], width=0.9, legend=False)
df_leakage.plot(ax=ax2, kind='bar', color=['gold', 'green', 'blue'], width=0.9, legend=False)
df_leakage.plot(ax=ax1, kind='bar', color=['gold', 'green', 'blue'], width=0.9)
for tick in ax3.get_xticklabels():
tick.set_rotation(90)
#d = .002
#kwargs = dict(transform=ax1.transAxes, color='k', clip_on=False)
#ax1.plot((-d, +d), (-d, +d), **kwargs)
#ax1.plot((1 - d, 1 + d), (-d, +d), **kwargs)
#kwargs.update(transform=ax2.transAxes)
#ax2.plot((-d, +d), (1 - d, 1 + d), **kwargs)
#ax2.plot((1 - d, 1 + d), (1 - d, 1 + d), **kwargs)
plt.subplots_adjust(hspace=0.15)
ax1.axvline(x=14.5, color='0.5', linewidth=2.0)
ax2.axvline(x=14.5, color='0.5', linewidth=2.0)
ax3.axvline(x=14.5, color='0.5', linewidth=2.0)
ax2.axhline(y=-0.5, color='0.5', linewidth=0.5, linestyle='--')
ax2.axhline(y=0.5, color='0.5', linewidth=0.5, linestyle='--')
#plt.title('energy leakage between netTOA and OHC')
ax2.set_ylabel('$W \; m^{-2}$')
plt.savefig('eei_leakage_breakdown.png', bbox_inches='tight', dpi=400)
plt.show()
In [ ]: