diff --git a/landice/output_processing_li/plot_mass_balance.py b/landice/output_processing_li/plot_mass_balance.py
index c7ee81773..7e6dcc582 100755
--- a/landice/output_processing_li/plot_mass_balance.py
+++ b/landice/output_processing_li/plot_mass_balance.py
@@ -30,61 +30,167 @@
 dyr[1:] = yr[1:]-yr[:-1]
 
 # ---- Figure: mass change rate ---
-fig = plt.figure(1, figsize=(6, 10), facecolor='w')
+fig, ax = plt.subplots(3,3, figsize=(15,12))
+#=====================
+# 1. Total mass budget
+#=====================
+# get fields and plot stuff
+vol = f.variables['totalIceVolume'][:]
+vol = vol * rhoi / 1.0e12
+dvol = np.zeros(yr.shape)
+dvol[1:] = (vol[1:]-vol[:-1]) / dyr[1:]
+
+# --- time-step-wise mass change ---
+ax[0,0].plot(yr, dvol, 'k', linewidth=3, label="mass change")
+
+SMB = f.variables['totalSfcMassBal'][:] / 1.0e12
+ax[0,0].plot(yr, SMB, label="SMB")
+
+BMB = f.variables['totalBasalMassBal'][:] / 1.0e12
+ax[0,0].plot(yr, BMB, label="BMB")
+
+calv = -f.variables['totalCalvingFlux'][:] / 1.0e12
+ax[0,0].plot(yr, calv, label='calv')
+
+FMF = -f.variables['totalFaceMeltingFlux'][:] / 1.0e12
+ax[0,0].plot(yr, FMF, label='facemelt')
+
+# Grounding line flux and grounding line migration flux
+# do not come into the global mass budget
+tot = SMB+BMB+calv+FMF
+ax[0,0].plot(yr, tot, "--", label="total")
+
+ax[0,0].legend(loc='best', prop={'size': 6})
+ax[0,0].set_ylabel('Mass change (Gt)')
+# --- cumulative mass change ---
 
-ax = fig.add_subplot(2, 1, 1)
-plt.xlabel('Year')
-plt.ylabel('Grounded mass change (Gt/yr)')
-plt.grid()
+ax[1,0].plot(yr, vol - vol[0], 'k', linewidth=3, label='total mass change')
+ax[1,0].plot(yr, (SMB*dyr).cumsum(), label="SMB")
+ax[1,0].plot(yr, (BMB*dyr).cumsum(), label="BMB")
 
+ax[1,0].plot(yr, (calv*dyr).cumsum(), label='calving')
+ax[1,0].plot(yr, (FMF*dyr).cumsum(), label='facemelt')
+ax[1,0].plot(yr, (tot * dyr).cumsum(), "--", label='total budget')
 
+ax[1,0].legend(loc='best', prop={'size': 6})
+ax[1,0].set_ylabel('Cumulative mass change (Gt)')
+
+ax[2,0].semilogy(yr, np.abs( (tot - dvol ) / dvol ) )
+ax[2,0].set_ylabel('fractional error', color='tab:blue')
+absErrAxTot = ax[2,0].twinx()
+absErrAxTot.plot(yr, tot - dvol, color='tab:orange')
+#========================
+# 2. Grounded mass budget
+#========================
 # get fields and plot stuff
 volGround = f.variables['groundedIceVolume'][:]
 volGround = volGround * rhoi / 1.0e12
 dvolGround = np.zeros(yr.shape)
 dvolGround[1:] = (volGround[1:]-volGround[:-1]) / dyr[1:]
-ax.plot(yr, dvolGround, 'k', linewidth=3, label="grounded mass change")
+ax[0,1].plot(yr, dvolGround, 'k', linewidth=3, label="mass change")
 
+# --- time-step-wise mass change ---
 SMBg = f.variables['totalGroundedSfcMassBal'][:] / 1.0e12
-ax.plot(yr, SMBg, label="grounded SMB")
+ax[0,1].plot(yr, SMBg, label="SMB")
 
 BMBg = f.variables['totalGroundedBasalMassBal'][:] / 1.0e12
-ax.plot(yr, BMBg, label="grounded BMB")
+ax[0,1].plot(yr, BMBg, label="BMB")
+
+calvg = -f.variables['totalGroundedCalvingFlux'][:] / 1.0e12
+ax[0,1].plot(yr, calvg, label='calv')
 
-GLflux = f.variables['groundingLineFlux'][:]/1.0e12
-ax.plot(yr, -GLflux, label="GL flux")
+FMFg = -f.variables['totalGroundedFaceMeltingFlux'][:] / 1.0e12
+ax[0,1].plot(yr, FMFg, label='facemelt')
 
-GLMigrationflux = f.variables['groundingLineMigrationFlux'][:]/1.0e12
-ax.plot(yr, -GLMigrationflux, label="GL migration flux")
+GLflux = -f.variables['groundingLineFlux'][:]/1.0e12
+ax[0,1].plot(yr, GLflux, label="GL flux")
 
-ax.plot(yr, SMBg+BMBg, "--", label="SMB+BMB")
-ax.plot(yr, SMBg+BMBg-GLflux, "--", label="SMB+BMB+GLf")
-ax.plot(yr, SMBg+BMBg-GLflux-GLMigrationflux, "--", label="SMB+BMB+GLf+GLMf")
+GLMigrationflux = -f.variables['groundingLineMigrationFlux'][:]/1.0e12
+ax[0,1].plot(yr, GLMigrationflux, label="GL migration flux")
 
-plt.legend(loc='best', prop={'size': 6})
-plt.tight_layout()
+grndTot = SMBg+calvg+BMBg+FMFg+GLflux+GLMigrationflux
+ax[0,1].plot(yr, grndTot, "--", label="total")
+
+ax[0,1].legend(loc='best', prop={'size': 6})
 
 # --- Figure: cumulative mass change ---
-#fig2 = plt.figure(2, figsize=(6,5), facecolor='w')
-ax = fig.add_subplot(2, 1, 2)
-plt.xlabel('Year')
-plt.ylabel('Cumulative grounded mass change (Gt)')
-plt.grid()
-
-plt.plot(yr, volGround - volGround[0], 'k', linewidth=3, label='grounded mass')
-plt.plot(yr, (SMBg*dyr).cumsum(), label="grounded SMB")
-plt.plot(yr, (BMBg*dyr).cumsum(), label="grounded BMB")
-plt.plot(yr, -1.0*(GLflux*dyr).cumsum(), label="GL flux")
-plt.plot(yr, -1.0*(GLMigrationflux*dyr).cumsum(), label="GL Migration flux")
-plt.plot(yr, ( (SMBg+BMBg) * dyr).cumsum(), "--", label='SMB+BMB')
-plt.plot(yr, ( (SMBg+BMBg-GLflux) * dyr).cumsum(), "--", label='SMB+BMB+GLflux')
-plt.plot(yr, ( (SMBg+BMBg-GLflux-GLMigrationflux) * dyr).cumsum(), "--", label='SMB+BMB+GLf+GLMf')
-#plt.plot(yr, dareag.cumsum(), 'k', label="mass change due to change in grounded area")
-
-plt.legend(loc='best', prop={'size': 6})
-plt.tight_layout()
-
-print("Plotting complete.")
 
+ax[1,1].plot(yr, volGround - volGround[0], 'k', linewidth=3, label='total mass change')
+ax[1,1].plot(yr, (SMBg*dyr).cumsum(), label="SMB")
+ax[1,1].plot(yr, (BMBg*dyr).cumsum(), label="BMB")
+
+ax[1,1].plot(yr, (calvg*dyr).cumsum(), label='calving')
+ax[1,1].plot(yr, (FMFg*dyr).cumsum(), label='facemelt')
+ax[1,1].plot(yr, (GLflux*dyr).cumsum(), label="GL flux")
+ax[1,1].plot(yr, (GLMigrationflux*dyr).cumsum(), label="GL Migration flux")
+ax[1,1].plot(yr, ( (grndTot) * dyr).cumsum(), "--", label='total budget')
+
+ax[1,1].legend(loc='best', prop={'size': 6})
+ax[2,1].semilogy(yr, np.abs( (grndTot - dvolGround)  / dvolGround ) )
+
+absErrAxGrnd = ax[2,1].twinx()
+absErrAxGrnd.plot(yr, grndTot - dvolGround, color='tab:orange')
+#========================
+# 3. Floating mass budget
+#========================
+
+# get fields and plot stuff
+volFloat = f.variables['floatingIceVolume'][:]
+volFloat = volFloat * rhoi / 1.0e12
+dvolFloat = np.zeros(yr.shape)
+dvolFloat[1:] = (volFloat[1:]-volFloat[:-1]) / dyr[1:]
+ax[0,2].plot(yr, dvolFloat, 'k', linewidth=3, label="mass change")
+
+SMBf = f.variables['totalFloatingSfcMassBal'][:] / 1.0e12
+ax[0,2].plot(yr, SMBf, label="SMB")
+
+BMBf = f.variables['totalFloatingBasalMassBal'][:] / 1.0e12
+ax[0,2].plot(yr, BMBf, label="BMB")
+
+calvf = -f.variables['totalFloatingCalvingFlux'][:] / 1.0e12
+ax[0,2].plot(yr, calvf, label='calv')
+
+FMFf = -f.variables['totalFloatingFaceMeltingFlux'][:] / 1.0e12
+ax[0,2].plot(yr, FMFf, label='facemelt')
+
+ax[0,2].plot(yr, -GLflux, label="GL flux")
+
+ax[0,2].plot(yr, -GLMigrationflux, label="GL migration flux")
+
+fltTot = SMBf+BMBf+FMFf+calvf-GLflux-GLMigrationflux
+ax[0,2].plot(yr, fltTot, "--", label="total")
+
+ax[0,2].legend(loc='best', prop={'size': 6})
+
+# --- cumulative mass change ---
+ax[1,2].plot(yr, volFloat - volFloat[0], 'k', linewidth=3, label='total mass change')
+ax[1,2].plot(yr, (SMBf*dyr).cumsum(), label="SMB")
+ax[1,2].plot(yr, (BMBf*dyr).cumsum(), label="BMB")
+
+ax[1,2].plot(yr, (calvf*dyr).cumsum(), label='calving')
+ax[1,2].plot(yr, (FMFf*dyr).cumsum(), label='facemelt')
+ax[1,2].plot(yr, (-GLflux*dyr).cumsum(), label="GL flux")
+ax[1,2].plot(yr, (-GLMigrationflux*dyr).cumsum(), label="GL Migration flux")
+ax[1,2].plot(yr, (fltTot * dyr).cumsum(), "--", label='total budget')
+
+ax[1,2].legend(loc='best', prop={'size': 6})
+
+ax[2,2].semilogy(yr, np.abs( (fltTot - dvolFloat) / dvolFloat ) )
+
+absErrAxFlt = ax[2,2].twinx()
+absErrAxFlt.plot(yr, fltTot - dvolFloat, color='tab:orange')
+absErrAxFlt.set_ylabel('Absolute error (Gt; budget - true)', color='tab:orange')
+
+for plotAx in ax.ravel():
+    plotAx.grid('on')
+ax[2,0].set_xlabel('yr')
+ax[2,1].set_xlabel('yr')
+ax[2,2].set_xlabel('yr')
+
+ax[0,0].set_title('Total budget')
+ax[0,1].set_title('Grounded budget')
+ax[0,2].set_title('Floating budget')
+
+fig.subplots_adjust(wspace=0.5)
 plt.show()
 f.close()