diff --git a/aice_animation.py b/aice_animation.py
index 285ef77..7d332b4 100644
--- a/aice_animation.py
+++ b/aice_animation.py
@@ -14,16 +14,16 @@
 # Directory containing CICE output files
 directory = '/short/y99/kaa561/roms_spinup_newest/cice/'
 # Number of time indices in each file
-num_ts = [73]
+num_ts = [18, 270, 270, 270, 252]
 # File number to start with for the animation (1-based)
-start_file = 1
+start_file = 5
 # Degrees to radians conversion factor
 deg2rad = pi/180
 # Names of each month for making titles
 month_names = ['January', 'February', 'March', 'April', 'May', 'June', 'July', 'August', 'September', 'October', 'November', 'December']
 
 # Read grid from the first file
-id = Dataset(directory + 'iceh1.nc', 'r')
+id = Dataset(directory + 'iceh' + str(start_file) + '.nc', 'r')
 lon = id.variables['TLON'][:,:]
 lat = id.variables['TLAT'][:,:]
 id.close()
@@ -62,10 +62,10 @@ def animate(i):
     img = contourf(x, y, data, lev, cmap='jet', extend='both')
     axis('off')
     # Add a title containing the date
-    title(str(time.day) + ' ' + month_names[time.month-1] + ' ' + str(time.year), fontsize=30)
+    title(str(time.day) + ' ' + month_names[time.month-1], fontsize=30) # + ' ' + str(time.year), fontsize=30)
     return img
 
 # Animate once every time index from start_file to the last file
-anim = FuncAnimation(fig, func=animate, frames=sum(array(num_ts[start_file-1:])))
+anim = FuncAnimation(fig, func=animate, frames=range(179,252)) #sum(array(num_ts[start_file-1:])))
 # Save as an mp4 with one frame per second
 anim.save('aice.mp4', fps=1)
diff --git a/circumpolar_plot.py b/circumpolar_plot.py
index 89215dd..4d86263 100644
--- a/circumpolar_plot.py
+++ b/circumpolar_plot.py
@@ -137,7 +137,7 @@ def circumpolar_plot (file_path, var_name, tstep, depth_key, depth, depth_bounds
             colour_map = 'jet'
     else:
         # Determine bounds automatically
-        if var_name in ['u', 'v', 'ubar', 'vbar', 'm', 'zeta', 'shflux', 'ssflux', 'sustr', 'svstr', 'bustr', 'bvstr']:
+        if var_name in ['u', 'v', 'ubar', 'vbar', 'm', 'shflux', 'ssflux', 'sustr', 'svstr', 'bustr', 'bvstr']:
             # Center levels on 0 for certain variables, with a blue-to-red
             # colourmap
             max_val = amax(abs(data))
diff --git a/dpt_timeseries.py b/dpt_timeseries.py
index 78a793a..9aef2a1 100644
--- a/dpt_timeseries.py
+++ b/dpt_timeseries.py
@@ -91,6 +91,7 @@ def dpt_timeseries (file_path):
     plot(time, transport)
     xlabel('Years')
     ylabel('Drake Passage Transport (Sv)')
+    grid(True)
     show(block=False)
     #savefig('dpt_timeseries.png')
 
diff --git a/file_guide.txt b/file_guide.txt
index 0fe535a..9cdcc1d 100644
--- a/file_guide.txt
+++ b/file_guide.txt
@@ -108,18 +108,18 @@ plot_volume.py: Extracts the volume values written in the ocean.log files (you
 spinup_plots.py: Analyse a ROMS spinup by calculating and plotting 9 timeseries:
                  total heat content, total salt content, area-averaged ice shelf
 		 melt rate, ice shelf basal mass loss, total kinetic energy,
-		 maximum velocity, Drake Passage transport, total sea ice
-		 extent, and average bottom water temperature in ice shelf
-		 cavities. Also write the timeseries to a log file so they don't
-		 have to be recomputed if the run is extended; at the beginning
-		 of this script, previous values will be read from the same
-		 log file if it exists. Note that to run this file, density
-		 anomalies (variable name "rho") must be written to the ocean
-		 history/averages file. This option can be activated in the .in
-		 ROMS configuration file with Hout(idDano) and Aout(idDano).
-		 Otherwise, you can use make_density_file.py to calculate
-		 absolute density based on the seawater equation of state;
-		 however, this process is quite slow.
+		 Drake Passage transport, total sea ice extent, net sea ice-to-
+		 ocean freshwater flux, and area-averaged bottom water
+		 temperature in ice shelf cavities. Also write the timeseries to
+		 a log file so they don't have to be recomputed if the run is
+		 extended; at the beginning of this script, previous values will
+		 be read from the same log file if it exists. Note that to run
+		 this file, density anomalies (variable name "rho") must be
+		 written to the ocean history/averages file. This option can be
+		 activated in the .in ROMS configuration file with Hout(idDano)
+		 and Aout(idDano). Otherwise, you can use make_density_file.py
+		 to calculate absolute density based on the seawater equation of
+		 state; however, this process is quite slow.
 		 To run: First edit the i and j coordinates of the Drake Passage
 		         (in the function calc_drakepsgtrans) for your grid, and
 			 the value of rho0 (in the function get_rho) to match
diff --git a/h_circumpolar.py b/h_circumpolar.py
index b7b0d54..df9690d 100644
--- a/h_circumpolar.py
+++ b/h_circumpolar.py
@@ -32,7 +32,7 @@ def h_circumpolar (grid_path, fig_name):
     x = -(lat+90)*cos(lon*deg2rad+pi/2)
     y = (lat+90)*sin(lon*deg2rad+pi/2)
 
-    lev = arange(0,6000,500)
+    lev = arange(0,max(data),500)
 
     # Plot
     fig = figure(figsize=(16,12))
@@ -43,7 +43,8 @@ def h_circumpolar (grid_path, fig_name):
     title('Bathymetry (m)', fontsize=30)
     axis('off')
 
-    savefig(fig_name)
+    show()
+    #savefig(fig_name)
 
 
 # Command-line interface
diff --git a/massloss.py b/massloss.py
index a793460..0a0cf26 100644
--- a/massloss.py
+++ b/massloss.py
@@ -25,7 +25,13 @@ def massloss (file_path, log_path):
     i_max = [350, 872, 975, 1030, 1090, 1155, 1190, 1240, 1369, 1443, 12]
     j_min = [1,   20,  150, 140,  160,  150,  187,  1,    65,   80,   100]
     j_max = [125, 123, 175, 160,  185,  200,  220,  135,  116,  150,  120]
-      
+    # Observed mass loss (Rignot 2013) and uncertainty for each ice shelf
+    obs_massloss = [35.5, 47.7, 144.9, 101.2, 51.8, 89, 20.7, 155.4, 9.7, 22.5]
+    obs_massloss_error = [23, 34, 14, 8, 19, 17, 67, 45, 16, 12]
+    # Observed ice shelf melt rates and uncertainty
+    obs_ismr = [0.6, 0.1, 4.3, 16.2, 1.7, 3.8, 0.4, 0.3, 0.1, 0.5]
+    obs_ismr_error = [0.4, 0.1, 0.4, 1, 0.6, 0.7, 1, 0.1, 0.2, 0.2]
+          
     # Density of ice in kg/m^3
     rho_ice = 916
 
@@ -82,7 +88,7 @@ def massloss (file_path, log_path):
 
     # Set up array of conversion factors from mass loss to area-averaged melt
     # rate for each ice shelf
-    factors = empty(12)
+    factors = empty(len(names))
 
     # Process each timestep separately to prevent memory overflow
     for t in range(start_t, size(time)):
@@ -120,18 +126,42 @@ def massloss (file_path, log_path):
     # Plot each timeseries
     print 'Plotting'
     for index in range(len(names)):
+        # Calculate the bounds on observed mass loss and melt rate
+        massloss_low = obs_massloss[index] - obs_massloss_error[index]
+        massloss_high = obs_massloss[index] + obs_massloss_error[index]
+        ismr_low = obs_ismr[index] - obs_ismr_error[index]
+        ismr_high = obs_ismr[index] + obs_ismr_error[index]
+        # Set up plot: mass loss and melt rate are directly proportional (with
+        # a different constant of proportionality for each ice shelf depending
+        # on its area) so plot one line with two y-axes
         fig, ax1 = subplots()
-        # Mass loss and melt rate are directly proportional (with a different
-        # constant of proportionality for each ice shelf depending on its area)
-        # so plot one line with two y-axes scales.
-        ax1.plot(time, massloss[index,:])        
+        ax1.plot(time, massloss[index,:], color='black')
+        # In blue, add dashed lines for observed mass loss
+        ax1.axhline(massloss_low, color='b', linestyle='dashed')
+        ax1.axhline(massloss_high, color='b', linestyle='dashed')
+        # Make sure y-limits won't cut off observed melt rate
+        ymin = min(0.95*ismr_low/factors[index], ax1.get_ylim()[0])
+        ymax = max(1.05*ismr_high/factors[index], ax1.get_ylim()[1])
+        ax1.set_ylim([ymin, ymax])
+        # Title and ticks in blue for this side of the plot
+        ax1.set_ylabel('Basal Mass Loss (Gt/y)', color='b')
+        for t1 in ax1.get_yticklabels():
+            t1.set_color('b')
         ax1.set_xlabel('Years')
-        ax1.set_ylabel('Basal Mass Loss (Gt/y)')
+        ax1.grid(True)
+        # Twin axis for melt rates
         ax2 = ax1.twinx()
-        ax2.set_ylabel('Area-Averaged Ice Shelf Melt Rate (m/y)')
         # Make sure the scales line up
-        limits = ax1.get_ylim()        
+        limits = ax1.get_ylim()
         ax2.set_ylim([limits[0]*factors[index], limits[1]*factors[index]])
+        # In red, add dashed lines for observed ice shelf melt rates
+        ax2.axhline(ismr_low, color='r', linestyle='dashed')
+        ax2.axhline(ismr_high, color='r', linestyle='dashed')
+        # Title and ticks in red for this side of the plot
+        ax2.set_ylabel('Area-Averaged Ice Shelf Melt Rate (m/y)', color='r')
+        for t2 in ax2.get_yticklabels():
+            t2.set_color('r')
+        # Name of the ice shelf for the main title
         title(names[index])
         fig.savefig(fig_names[index])
         
diff --git a/plot_kinetic_energy.py b/plot_kinetic_energy.py
index 4f674b4..0259eec 100644
--- a/plot_kinetic_energy.py
+++ b/plot_kinetic_energy.py
@@ -69,6 +69,7 @@ def plot_kinetic_energy (files, dt, freq):
     plot(time, ke_all)
     xlabel('Years')
     ylabel('Kinetic energy')
+    grid(True)
     show()
 
 
diff --git a/plot_maxspeed.py b/plot_maxspeed.py
index 5ba5ca1..32deede 100644
--- a/plot_maxspeed.py
+++ b/plot_maxspeed.py
@@ -75,6 +75,7 @@ def plot_maxspeed (files, dt, freq):
     plot(time, maxv_all)
     xlabel('Years')
     ylabel('Maximum speed')
+    grid(True)
     show(block=False)
 
 
diff --git a/plot_volume.py b/plot_volume.py
index d639517..d01c14f 100644
--- a/plot_volume.py
+++ b/plot_volume.py
@@ -73,6 +73,7 @@ def plot_volume (files, dt, freq):
     plot(time, vol_all)
     xlabel('Years')
     ylabel('Volume Percent Anomalies')
+    grid(True)
     show()
 
 
diff --git a/repeat_forcing.py b/repeat_forcing.py
index 416567a..9fa8921 100644
--- a/repeat_forcing.py
+++ b/repeat_forcing.py
@@ -74,12 +74,12 @@ def process (directory, head, tail, year_start, perday, var_list):
 # User parameters to edit here
 
 # Data where files <an_head>yyyy<tail> and <fc_head>yyyy<tail> exist
-directory = '../data/ERA_Interim/'
+directory = '/short/m68/kaa561/ROMS-CICE-MCT/data/ERA_Interim/subdaily/'
 # First part of filename for AN and FC files
 an_head = 'AN_'
 fc_head = 'FC_'
 # Last part of filename (common to both AN and FC)
-tail = '_unlim.nc'
+tail = '_subdaily.nc'
 # Year to build data from
 year_start = 1995
 # Number of records per day
diff --git a/romscice_atm_subdaily.py b/romscice_atm_subdaily.py
index 44708a4..22d9799 100644
--- a/romscice_atm_subdaily.py
+++ b/romscice_atm_subdaily.py
@@ -3,17 +3,19 @@
 from scipy.interpolate import LinearNDInterpolator, RectBivariateSpline
 
 # Convert two ERA-Interim files:
-# AN_yyyy_unlim_orig.nc: one year of 6-hour measurements for surface pressure 
-#                        (sp), 2-metre temperature (t2m) and dew point (d2m),
-#                        total cloud cover (tcc), and 10-metre winds (u10, v10)
-# FC_yyyy_unlim_orig.nc: one year of 12-hour measurements for total 
-#                        precipitation (tp) 
+# AN_yyyy_subdaily_orig.nc: one year of 6-hour measurements for surface pressure
+#                           (sp), 2-metre temperature (t2m) and dew point (d2m),
+#                           total cloud cover (tcc), and 10-metre winds (u10, 
+#                           v10)
+# FC_yyyy_subdaily_orig.nc: one year of 12-hour measurements for total 
+#                           precipitation (tp) and snowfall (sf) 
 # to two ROMS-CICE input forcing files with the correct units and naming 
 # conventions:
-# AN_yyyy_unlim.nc: one year of 6-hour measurements for surface pressure
+# AN_yyyy_subdaily.nc: one year of 6-hour measurements for surface pressure
 #                   (Pair), temperature (Tair), specific humidity (Qair), 
 #                   cloud fraction (cloud), and winds (Uwind, Vwind)
-# FC_yyyy_unlim.nc: one year of 12-hour measurements for rainfall (rain)
+# FC_yyyy_subdaily.nc: one year of 12-hour measurements for rainfall (rain) and
+#                   snowfall (snow)
 # Input: year = integer containing the year to process
 #        count = time record in the given year to start with
 
@@ -30,11 +32,11 @@ def convert_file (year, count):
 
     # Paths of ROMS grid file, input ERA-Interim files, and output ROMS-CICE
     # files; other users will need to change these
-    grid_file = '/short/m68/kaa561/roms_circumpolar/data/caisom001_OneQuartergrd.nc'
-    input_atm_file = '../data/ERA_Interim/AN_' + str(year) + '_unlim_orig.nc'
-    input_ppt_file = '../data/ERA_Interim/FC_' + str(year) + '_unlim_orig.nc'
-    output_atm_file = '../data/ERA_Interim/AN_' + str(year) + '_unlim.nc'
-    output_ppt_file = '../data/ERA_Interim/FC_' + str(year) + '_unlim.nc'
+    grid_file = '/short/m68/kaa561/ROMS-CICE-MCT/apps/common/grid/circ30S_quarterdegree_rp5.nc'
+    input_atm_file = '/short/m68/kaa561/ROMS-CICE-MCT/data/ERA_Interim/subdaily/AN_' + str(year) + '_subdaily_orig.nc'
+    input_ppt_file = '/short/m68/kaa561/ROMS-CICE-MCT/data/ERA_Interim/subdaily/FC_' + str(year) + '_subdaily_orig.nc'
+    output_atm_file = '/short/m68/kaa561/ROMS-CICE-MCT/data/ERA_Interim/subdaily/AN_' + str(year) + '_subdaily.nc'
+    output_ppt_file = '/short/m68/kaa561/ROMS-CICE-MCT/data/ERA_Interim/subdaily/FC_' + str(year) + '_subdaily.nc'
     logfile = str(year) + '.log'
 
     Lv = 2.5e6 # Latent heat of vapourisation, J/kg
@@ -139,6 +141,9 @@ def convert_file (year, count):
         tair = interp_era2roms(t2m, lon_era, lat_era, lon_roms, lat_roms)
         oatm_fid.variables['Tair'][t,:,:] = tair-273.15
         qair = interp_era2roms(rh, lon_era, lat_era, lon_roms, lat_roms)
+        # Constrain humidity values to be between 0 and 1
+        qair[qair < 0] = 0.0
+        qair[qair > 1] = 1.0
         oatm_fid.variables['Qair'][t,:,:] = qair
         cloud = interp_era2roms(tcc, lon_era, lat_era, lon_roms, lat_roms)
         # Constrain cloud fractions to be between 0 and 1
@@ -183,6 +188,9 @@ def convert_file (year, count):
         oppt_fid.createVariable('rain', 'f8', ('time', 'eta_rho', 'xi_rho'))
         oppt_fid.variables['rain'].long_name = 'rain fall rate'
         oppt_fid.variables['rain'].units = 'm_per_12hr'
+        oppt_fid.createVariable('snow', 'f8', ('time', 'eta_rho', 'xi_rho'))
+        oppt_fid.variables['snow'].long_name = 'snow fall rate'
+        oppt_fid.variables['snow'].units = 'm_per_12hr'
         oppt_fid.close()
 
     log = open(logfile, 'a')
@@ -198,15 +206,19 @@ def convert_file (year, count):
         log.close()
         # Write the current time value to output FC file
         oppt_fid.variables['time'][t] = ppt_time[t]
-        # Read rainfall for this timestep
+        # Read data for this timestep
         ippt_fid = Dataset(input_ppt_file, 'r')
         tp = transpose(ippt_fid.variables['tp'][t,:,:])
+        sf = transpose(ippt_fid.variables['sf'][t,:,:])
         ippt_fid.close()
         # Interpolate to ROMS grid and write to output FC file
         rain = interp_era2roms(tp, lon_era, lat_era, lon_roms, lat_roms)
+        snow = interp_era2roms(sf, lon_era, lat_era, lon_roms, lat_roms)
         # Make sure there are no negative values
         rain[rain < 0] = 0.0
         oppt_fid.variables['rain'][t,:,:] = rain
+        snow[snow < 0] = 0.0
+        oppt_fid.variables['snow'][t,:,:] = snow
         oppt_fid.close()
 
     log = open(logfile, 'a')
diff --git a/spinup_plots.py b/spinup_plots.py
index b25a341..e8741be 100644
--- a/spinup_plots.py
+++ b/spinup_plots.py
@@ -4,13 +4,12 @@
 from os.path import *
 from cartesian_grid_3d import *
 
-# Analyse a ROMS spinup by calculating and plotting 10 timeseries:
+# Analyse a ROMS spinup by calculating and plotting 9 timeseries:
 # Total heat content
 # Total salt content
 # Area-averaged ice shelf melt rate
 # Ice shelf basal mass loss
 # Total kinetic energy
-# Maximum velocity
 # Drake Passage transport
 # Total sea ice extent
 # Net sea ice-to-ocean freshwater flux
@@ -126,33 +125,13 @@ def calc_totalsalt (file_path, dV, rho, t):
     return totalsalt    
 
 
-# Calculate area-averaged ice shelf melt rate at the given timestep t.
-# Input:
-# file_path = path to ocean history/averages file
-# dA = elements of area on the rho grid, masked with zice
-# t = timestep index in file_path
-# Output: avgismr = area-averaged ice shelf melt rate (m/y)
-#         ismr = 2D ice shelf melt rate field (m/y) at this timestep
-def calc_avgismr (file_path, dA, t):
-
-    # Read ice shelf melt rate, converting to float128 to prevent overflow
-    # during integration
-    id = Dataset(file_path, 'r')
-    ismr = ma.asarray(id.variables['m'][t,:,:], dtype=float128)
-    # Convert from m/s to m/y
-    ismr = ismr*365.25*24*60*60
-    id.close()    
-
-    # Integrate ismr over area and divide by total area to get average
-    avgismr = sum(ismr*dA)/sum(dA)
-    return avgismr, ismr
-
-
 # Calculate net basal mass loss based on the given ice shelf melt rate field.
 # Input:
 # ismr = 2D ice shelf melt rate field (m/y)
 # dA = elements of area on the rho grid, masked with zice
-# Output: massloss = net basal mass loss (Gt/y)
+# Output: 
+# massloss = net basal mass loss (Gt/y)
+# massloss_factor = conversion factor from mass loss to area-averaged melt rate
 def calc_massloss (ismr, dA):
 
     # Density of ice in kg/m^3
@@ -162,7 +141,8 @@ def calc_massloss (ismr, dA):
     volumeloss = sum(ismr*dA)
     # Convert to mass loss in Gt/y
     massloss = 1e-12*rho_ice*volumeloss
-    return massloss
+    massloss_factor = 1e12/(rho_ice*sum(dA))
+    return massloss, massloss_factor
 
 
 # Calculate total kinetic energy at the given timestep t.
@@ -195,15 +175,7 @@ def calc_tke (file_path, dV, rho, t):
 
     # Integrate 0.5*rho*(u^2 + v^2) over volume to get TKE
     tke = sum(0.5*rho*(u_rho**2 + v_rho**2)*dV)
-    return tke, u_rho, v_rho
-
-
-# Calculate the maximum velocity.
-# Input: u_rho, v_rho = u and v at timestep t, interpolated to the rho-grid
-# Output: maxvel = maximum velocity (m/s)
-def calc_maxvel (u_rho, v_rho):
-
-    return amax(sqrt(u_rho**2 + v_rho**2))
+    return tke
 
 
 # Calculate zonal transport through the Drake Passage.
@@ -336,13 +308,18 @@ def calc_bwtemp (file_path, dA, t):
 #            calculated values following computation)
 def spinup_plots (file_path, cice_path, log_path):
 
+    # Observed basal mass loss (Rignot 2013) and uncertainty
+    obs_massloss = 1325
+    obs_massloss_error = 235
+    # Observed ice shelf melt rates and uncertainty
+    obs_ismr = 0.85
+    obs_ismr_error = 0.1
+
     time = []
     ohc = []
     totalsalt = []
-    avgismr = []
     massloss = []
     tke = []
-    maxvel = []
     drakepsgtrans = []
     totalice = []
     totalfwflux = []
@@ -369,11 +346,6 @@ def spinup_plots (file_path, cice_path, log_path):
                 totalsalt.append(float(line))
             except (ValueError):
                 break
-        for line in f:
-            try:
-                avgismr.append(float(line))
-            except (ValueError):
-                break
         for line in f:
             try:
                 massloss.append(float(line))
@@ -384,11 +356,6 @@ def spinup_plots (file_path, cice_path, log_path):
                 tke.append(float(line))
             except (ValueError):
                 break
-        for line in f:
-            try:
-                maxvel.append(float(line))
-            except (ValueError):
-                break
         for line in f:
             try:
                 drakepsgtrans.append(float(line))
@@ -427,16 +394,11 @@ def spinup_plots (file_path, cice_path, log_path):
         ohc.append(calc_ohc(file_path, dV, rho, t))
         print 'Calculating total salt content'
         totalsalt.append(calc_totalsalt(file_path, dV, rho, t))
-        print 'Calculating average ice shelf melt rate'
-        avgismr_tmp, ismr = calc_avgismr(file_path, dA, t)
-        avgismr.append(avgismr_tmp)
         print 'Calculating basal mass loss'
-        massloss.append(calc_massloss(ismr, dA))
+        massloss_tmp, massloss_factor = calc_massloss(ismr, dA)
+        massloss.append(massloss_tmp)
         print 'Calculating total kinetic energy'
-        tke_tmp, u_rho, v_rho = calc_tke(file_path, dV, rho, t)
-        tke.append(tke_tmp)
-        print 'Calculating maximum velocity'
-        maxvel.append(calc_maxvel(u_rho, v_rho))
+        tke.append(calc_tke(file_path, dV, rho, t))
         print 'Calculating Drake Passage transport'
         drakepsgtrans.append(calc_drakepsgtrans(file_path, dy_wct, t))
         print 'Calculating total sea ice extent'
@@ -452,49 +414,79 @@ def spinup_plots (file_path, cice_path, log_path):
     plot(time, ohc)
     xlabel('Years')
     ylabel('Southern Ocean Heat Content (J)')
+    grid(True)
     savefig('ohc.png')
+
     print 'Plotting total salt content'
     clf()
     plot(time, totalsalt)
     xlabel('Years')
     ylabel('Southern Ocean Salt Content (kg)')
+    grid(True)
     savefig('totalsalt.png')
-    print 'Plotting average ice shelf melt rate'
-    clf()
-    plot(time, avgismr)
-    xlabel('Years')
-    ylabel('Area-averaged Ice Shelf Melt Rate (m/y)')
-    savefig('avgismr.png')
-    print 'Plotting basal mass loss'
+
+    print 'Plotting basal mass loss and area-averaged ice shelf melt rate'
     clf()
-    plot(time, massloss)
-    xlabel('Years')
-    ylabel('Ice Shelf Basal Mass Loss (Gt/y)')
+    # Calculate the bounds on observed mass loss and melt rate
+    massloss_low = obs_massloss - obs_massloss_error
+    massloss_high = obs_massloss + obs_massloss_error
+    ismr_low = obs_ismr - obs_ismr_error
+    ismr_high = obs_ismr + obs_ismr_error
+    # Set up plot: mass loss and melt rate are directly proportional, so plot
+    # one line with two y-axes
+    ax1 = subplot(111)
+    ax1.plot(time, massloss, color='black')
+    # In blue, add dashed lines for observed mass loss
+    ax1.axhline(massloss_low, color='b', linestyle='dashed')
+    ax1.axhline(massloss_high, color='b', linestyle='dashed')
+    # Make sure y-limits won't cut off observed melt rate
+    ymin = min(0.95*ismr_low/massloss_factor, ax1.get_ylim()[0])
+    ymax = max(1.05*ismr_high/massloss_factor, ax1.get_ylim()[1])
+    ax1.set_ylim([ymin, ymax])
+    # Title and ticks in blue for this side of the plot
+    ax1.set_ylabel('Basal Mass Loss (Gt/y)', color='b')
+    for t1 in ax1.get_yticklabels():
+        t1.set_color('b')        
+    ax1.set_xlabel('Years')
+    ax1.grid(True)
+    # Twin axis for melt rates
+    ax2 = ax1.twinx()
+    # Make sure the scales line up
+    limits = ax1.get_ylim()
+    ax2.set_ylim([limits[0]*massloss_factor, limits[1]*massloss_factor])
+    # In red, add dashed lines for observed ice shelf melt rates
+    ax2.axhline(ismr_low, color='r', linestyle='dashed')
+    ax2.axhline(ismr_high, color='r', linestyle='dashed')
+    # Title and ticks in red for this side of the plot
+    ax2.set_ylabel('Area-averaged Ice Shelf Melt Rate (m/y)', color='r')
+    for t2 in ax2.get_yticklabels():
+        t2.set_color('r')
     savefig('massloss.png')
+
     print 'Plotting total kinetic energy'
     clf()
     plot(time, tke)
     xlabel('Years')
     ylabel('Southern Ocean Total Kinetic Energy (J)')
+    grid(True)
     savefig('tke.png')
-    print 'Plotting maximum velocity'
-    clf()
-    plot(time, maxvel)
-    xlabel('Years')
-    ylabel('Maximum Southern Ocean Velocity (m/s)')
-    savefig('maxvel.png')
+
     print 'Plotting Drake Passage transport'
     clf()
     plot(time, drakepsgtrans)
     xlabel('Years')
     ylabel('Drake Passage Transport (Sv)')
+    grid(True)
     savefig('drakepsgtrans.png')
+
     print 'Plotting total sea ice extent'
     clf()
     plot(time, totalice)
     xlabel('Years')
     ylabel(r'Total Sea Ice Extent (million km$^2$)')
+    grid(True)
     savefig('totalice.png')
+
     print 'Plotting total sea ice-to-ocean freshwater flux'
     clf()
     plot(time, totalfwflux)
@@ -502,12 +494,15 @@ def spinup_plots (file_path, cice_path, log_path):
     plot(time, zeros(len(totalfwflux)), color='black')
     xlabel('Years')
     ylabel('Sea Ice-to-Ocean Freshwater Flux (Sv)')
+    grid(True)
     savefig('totalfwflux.png')
+
     print 'Plotting bottom water temperature'
     clf()
     plot(time, bwtemp)
     xlabel('Years')
     ylabel(r'Average Bottom Water Temperature in Ice Shelf Cavities ($^{\circ}$C)')
+    grid(True)
     savefig('bwtemp.png')
 
     print 'Saving results to log file'
@@ -521,18 +516,12 @@ def spinup_plots (file_path, cice_path, log_path):
     f.write('Southern Ocean Salt Content (kg):\n')
     for elm in totalsalt:
         f.write(str(elm) + '\n')
-    f.write('Area-averaged Ice Shelf Melt Rate (m/y):\n')
-    for elm in avgismr:
-        f.write(str(elm) + '\n')
     f.write('Ice Shelf Basal Mass Loss (Gt/y):\n')
     for elm in massloss:
         f.write(str(elm) + '\n')
     f.write('Southern Ocean Total Kinetic Energy (J):\n')
     for elm in tke:
         f.write(str(elm) + '\n')
-    f.write('Maximum Southern Ocean Velocity (m/s):\n')
-    for elm in maxvel:
-        f.write(str(elm) + '\n')
     f.write('Drake Passage Transport (Sv):\n')
     for elm in drakepsgtrans:
         f.write(str(elm) + '\n')