diff --git a/bugs_acc_fig.py b/bugs_acc_fig.py
new file mode 100644
index 0000000..94cd783
--- /dev/null
+++ b/bugs_acc_fig.py
@@ -0,0 +1,82 @@
+from netCDF4 import Dataset, num2date
+from numpy import *
+from matplotlib.pyplot import *
+from rotate_vector_roms import *
+
+def bugs_acc_fig (grid_path, laplacian_file, biharmonic_file, tstep):
+
+    # Month names for title
+    month_names = ['January', 'February', 'March', 'April', 'May', 'June', 'July', 'August', 'September', 'October', 'November', 'December']
+    # Degrees to radians conversion factor
+    deg2rad = pi/180.0
+
+    # Read angle and grid
+    id = Dataset(grid_path, 'r')
+    angle = id.variables['angle'][:-15,:]
+    lon = id.variables['lon_rho'][:-15,:-1]
+    lat = id.variables['lat_rho'][:-15,:-1]
+    id.close()
+
+    # Read surface velocity
+    id = Dataset(laplacian_file, 'r')
+    ur_lap = id.variables['u'][tstep-1,-1,:-15,:]
+    vr_lap = id.variables['v'][tstep-1,-1,:-15,:]
+    # Rotate to lon-lat space
+    u_lap, v_lap = rotate_vector_roms(ur_lap, vr_lap, angle)
+    # Now that they're on the same grid, get the speed
+    speed_laplacian = sqrt(u_lap**2 + v_lap**2)    
+    # Also read time and convert to Date object
+    time_id = id.variables['ocean_time']
+    time = num2date(time_id[tstep-1], units=time_id.units, calendar=time_id.calendar.lower())
+    # Get the date for the title
+    date_string = str(time.day) + ' ' + month_names[time.month-1] + ' ' + str(time.year)
+    id.close()
+    # Repeat velocity for biharmonic simulation
+    id = Dataset(biharmonic_file, 'r')
+    ur_bih = id.variables['u'][tstep-1,-1,:-15,:]
+    vr_bih = id.variables['v'][tstep-1,-1,:-15,:]
+    u_bih, v_bih = rotate_vector_roms(ur_bih, vr_bih, angle)
+    speed_biharmonic = sqrt(u_bih**2 + v_bih**2)
+    id.close()
+
+    # Calculate x and y coordinates for plotting circumpolar projection
+    x = -(lat+90)*cos(lon*deg2rad+pi/2)
+    y  = (lat+90)*sin(lon*deg2rad+pi/2)
+
+    # Make the plot
+    fig = figure(figsize=(20,10))
+    # Laplacian
+    ax = fig.add_subplot(1,2,1, aspect='equal')
+    pcolor(x, y, speed_laplacian, vmin=0, vmax=1.25, cmap='cool')
+    title('Laplacian viscosity', fontsize=24)
+    xlim([-40, 0])
+    ylim([0, 40])
+    ax.set_xticks([])
+    ax.set_yticks([])
+    # Biharmonic
+    ax = fig.add_subplot(1,2,2, aspect='equal')
+    img = pcolor(x, y, speed_biharmonic, vmin=0, vmax=1.25, cmap='cool')
+    title('Biharmonic viscosity', fontsize=24)
+    xlim([-40, 0])
+    ylim([0, 40])
+    ax.set_xticks([])
+    ax.set_yticks([])
+    # Colourbar on the right
+    cbaxes = fig.add_axes([0.94, 0.3, 0.02, 0.4])
+    cbar = colorbar(img, cax=cbaxes, extend='max', ticks=arange(0,1.25+0.25,0.25))
+    cbar.ax.tick_params(labelsize=16)
+    suptitle('Surface speed (m/s): snapshot on ' + date_string, fontsize=30)
+    subplots_adjust(wspace=0.05)
+    fig.show()
+    fig.savefig('bugs_acc.png')
+
+
+# Command-line interface
+if __name__ == "__main__":
+
+    grid_path = raw_input("Path to ROMS grid file: ")
+    laplacian_file = raw_input("Path to ocean history file for Laplacian simulation: ")
+    biharmonic_file = raw_input("Path to same ocean history file for baseline simulation: ")
+    tstep = int(raw_input("Timestep to plot (starting at 1): "))
+    bugs_acc_fig(grid_path, laplacian_file, biharmonic_file, tstep)
+    
diff --git a/bugs_dpt_fig.py b/bugs_dpt_fig.py
new file mode 100644
index 0000000..4ec9993
--- /dev/null
+++ b/bugs_dpt_fig.py
@@ -0,0 +1,74 @@
+from netCDF4 import Dataset
+from numpy import *
+from matplotlib.pyplot import *
+from mip_timeseries import roms_annual_avg
+
+def bugs_dpt_fig (dpt_log_laplacian, dpt_log_biharmonic):
+
+    year_start = 1992
+    
+    # Read Laplacian timeseries
+    time = []
+    dpt_laplacian = []
+    f = open(dpt_log_laplacian, 'r')
+    # Skip first line (header for time array)
+    f.readline()
+    for line in f:
+        try:
+            time.append(float(line))
+        except(ValueError):
+            # Reached the header for the next variable
+            break
+    for line in f:
+        dpt_laplacian.append(float(line))
+    f.close()
+    # Read biharmonic timeseries
+    dpt_biharmonic = []
+    f = open(dpt_log_biharmonic, 'r')
+    f.readline()
+    # Skip the time array, as it's the same
+    for line in f:
+        try:
+            tmp = float(line)
+        except(ValueError):
+            break
+    for line in f:
+        dpt_biharmonic.append(float(line))
+    f.close()
+    # Make sure they're the same length
+    if len(dpt_laplacian) != len(dpt_biharmonic):
+        print 'These logfiles do not cover the same time periods!'
+        exit
+    # Add start year to time array
+    time = array(time) + year_start
+    # Annually average DPT
+    dpt_laplacian_avg = roms_annual_avg(time, dpt_laplacian, year_start)
+    dpt_biharmonic_avg = roms_annual_avg(time, dpt_biharmonic, year_start)
+    # Make new time array
+    time_annual = arange(len(dpt_laplacian_avg)) + year_start
+    # Plot
+    fig, ax = subplots(figsize=(10,6))
+    ax.plot(time_annual, dpt_laplacian_avg, label='Laplacian', color='green', linewidth=2)
+    ax.plot(time_annual, dpt_biharmonic_avg, label='Biharmonic', color='blue', linewidth=2)
+    title('Drake Passage Transport (annually averaged)', fontsize=18)
+    xlabel('Year', fontsize=14)
+    ylabel('Sv', fontsize=14)
+    xlim([time_annual[0], time_annual[-1]])
+    grid(True)
+    # Move plot over to make room for legend
+    box = ax.get_position()
+    ax.set_position([box.x0, box.y0, box.width*0.8, box.height])
+    # Make legend
+    ax.legend(loc='center left', bbox_to_anchor=(1,0.5))
+    fig.show()
+    fig.savefig('bugs_dpt.png')
+
+
+# Command-line interface
+if __name__ == "__main__":
+
+    dpt_log_laplacian = raw_input("Path to Drake Passage Transport logfile for Laplacian viscosity simulation: ")
+    dpt_log_biharmonic = raw_input("Path to Drake Passage Transport logfile for baseline simulation: ")
+    bugs_dpt_fig(dpt_log_laplacian, dpt_log_biharmonic)
+
+    
diff --git a/bugs_drift_slices.py b/bugs_drift_slices.py
new file mode 100644
index 0000000..237df5d
--- /dev/null
+++ b/bugs_drift_slices.py
@@ -0,0 +1,228 @@
+from netCDF4 import Dataset
+from numpy import *
+from matplotlib.pyplot import *
+from calc_z import *
+from interp_lon_roms import *
+
+def bugs_drift_slices (grid_file, upwind_file, akima_file, split_file):
+
+    # Paths to ECCO2 files with initial conditions for temp and salt
+    ecco_temp_file = '/short/m68/kaa561/metroms_iceshelf/data/originals/ECCO2/THETA.1440x720x50.199201.nc'
+    ecco_salt_file = '/short/m68/kaa561/metroms_iceshelf/data/originals/ECCO2/SALT.1440x720x50.199201.nc'
+    # Longitude to interpolate to (OE)
+    lon0 = 0
+    # Bounds on plot
+    lat_min = -73
+    lat_max = -30
+    depth_min = -6000
+    depth_max = 0
+    # ROMS grid parameters
+    theta_s = 7.0
+    theta_b = 2.0
+    hc = 250
+    N = 31
+    # Bounds on colour scales for temperature and salinity
+    temp_min = -2
+    temp_max = 6
+    salt_min = 33.9
+    salt_max = 34.9
+    # Contours to overlay
+    temp_contour = 0.75
+    salt_contour = 34.5
+
+    # Get longitude for the title
+    if lon0 < 0:
+        lon_string = str(int(round(-lon0))) + r'$^{\circ}$W'
+    else:
+        lon_string = str(int(round(lon0))) + r'$^{\circ}$E'
+
+    print 'Processing ECCO2'
+    id = Dataset(ecco_temp_file, 'r')
+    # Read grid variables
+    ecco_lat = id.variables['LATITUDE_T'][:]
+    ecco_depth = -1*id.variables['DEPTH_T'][:]
+    if lon0 == 0:
+        # Hard-coded lon0 = 0E: average between the first (0.125 E) and last
+        # (359.875 E = -0.125 W) indices in the regular ECCO2 grid
+        ecco_temp = 0.5*(id.variables['THETA'][0,:,:,0] + id.variables['THETA'][0,:,:,-1])
+        id.close()
+        id = Dataset(ecco_salt_file, 'r')
+        ecco_salt = 0.5*(id.variables['SALT'][0,:,:,0] + id.variables['SALT'][0,:,:,-1])
+        id.close()
+    else:
+        print 'lon0 is only coded for 0E at this time'
+        return
+
+    print 'Building ROMS grid'
+    id = Dataset(grid_file, 'r')
+    lon_2d = id.variables['lon_rho'][:,:]
+    lat_2d = id.variables['lat_rho'][:,:]
+    h = id.variables['h'][:,:]
+    zice = id.variables['zice'][:,:]
+    id.close()
+    # Get a 3D array of z-coordinates; sc_r and Cs_r are unused in this script
+    z_3d, sc_r, Cs_r = calc_z(h, zice, theta_s, theta_b, hc, N)
+    # Make sure we are in the range 0-360
+    if lon0 < 0:
+        lon0 += 360
+    print 'Processing upwind advection'
+    id = Dataset(upwind_file, 'r')
+    upwind_temp_3d = id.variables['temp'][0,:,:,:]
+    upwind_salt_3d = id.variables['salt'][0,:,:,:]
+    id.close()
+    # Interpolate to lon0
+    upwind_temp, z, lat = interp_lon_roms(upwind_temp_3d, z_3d, lat_2d, lon_2d, lon0)
+    upwind_salt, z, lat = interp_lon_roms(upwind_salt_3d, z_3d, lat_2d, lon_2d, lon0)
+    print 'Processing Akima advection'
+    id = Dataset(akima_file, 'r')
+    akima_temp_3d = id.variables['temp'][0,:,:,:]
+    akima_salt_3d = id.variables['salt'][0,:,:,:]
+    id.close()
+    akima_temp, z, lat = interp_lon_roms(akima_temp_3d, z_3d, lat_2d, lon_2d, lon0)
+    akima_salt, z, lat = interp_lon_roms(akima_salt_3d, z_3d, lat_2d, lon_2d, lon0)    
+    print 'Processing split advection'
+    id = Dataset(split_file, 'r')
+    split_temp_3d = id.variables['temp'][0,:,:,:]
+    split_salt_3d = id.variables['salt'][0,:,:,:]
+    id.close()
+    split_temp, z, lat = interp_lon_roms(split_temp_3d, z_3d, lat_2d, lon_2d, lon0)
+    split_salt, z, lat = interp_lon_roms(split_salt_3d, z_3d, lat_2d, lon_2d, lon0)
+    # Switch back to range -180-180
+    if lon0 > 180:
+        lon0 -= 360
+
+    # Set up axis labels the way we want them
+    lat_ticks = arange(lat_min+3, lat_max+10, 10)
+    lat_labels = []
+    for val in lat_ticks:
+        lat_labels.append(str(int(round(-val))) + r'$^{\circ}$S')
+    depth_ticks = range(depth_min+1000, 0+1000, 1000)
+    depth_labels = []
+    for val in depth_ticks:
+        depth_labels.append(str(int(round(-val))))
+
+    print 'Plotting'
+    fig = figure(figsize=(14,24))
+    # ECCO2
+    gs1 = GridSpec(1,2)
+    gs1.update(left=0.1, right=0.95, bottom=0.7575, top=0.94, wspace=0.08)
+    # Temperature
+    ax = subplot(gs1[0,0])
+    pcolor(ecco_lat, ecco_depth, ecco_temp, vmin=temp_min, vmax=temp_max, cmap='jet')
+    # Overlay contour
+    contour(ecco_lat, ecco_depth, ecco_temp, levels=[temp_contour], color='black')
+    title(r'Temperature ($^{\circ}$C)', fontsize=24)
+    ylabel('Depth (m)', fontsize=18)
+    xlim([lat_min, lat_max])
+    ylim([depth_min, depth_max])
+    ax.set_xticks(lat_ticks)
+    ax.set_xticklabels(lat_labels, fontsize=16)
+    ax.set_yticks(depth_ticks)
+    ax.set_yticklabels(depth_labels, fontsize=16)
+    text(-64, 1000, 'a) ECCO2 initial conditions at ' + lon_string + ', January 1992', fontsize=28)
+    # Salinity
+    ax = subplot(gs1[0,1])
+    pcolor(ecco_lat, ecco_depth, ecco_salt, vmin=salt_min, vmax=salt_max, cmap='jet')
+    contour(ecco_lat, ecco_depth, ecco_salt, levels=[salt_contour], color='black')
+    title('Salinity (psu)', fontsize=24)
+    xlim([lat_min, lat_max])
+    ylim([depth_min, depth_max])
+    ax.set_xticks(lat_ticks)
+    ax.set_xticklabels(lat_labels, fontsize=16)
+    ax.set_yticks(depth_ticks)
+    ax.set_yticklabels([])
+    # Upwind advection
+    gs2 = GridSpec(1,2)
+    gs2.update(left=0.1, right=0.95, bottom=0.525, top=0.7075, wspace=0.08)
+    # Temperature
+    ax = subplot(gs2[0,0])
+    pcolor(lat, z, upwind_temp, vmin=temp_min, vmax=temp_max, cmap='jet')
+    contour(lat, z, upwind_temp, levels=[temp_contour], color='black')
+    ylabel('Depth (m)', fontsize=18)
+    xlim([lat_min, lat_max])
+    ylim([depth_min, depth_max])
+    ax.set_xticks(lat_ticks)
+    ax.set_xticklabels(lat_labels, fontsize=16)
+    ax.set_yticks(depth_ticks)
+    ax.set_yticklabels(depth_labels, fontsize=16)
+    text(-60, 300, 'b) Upwind third-order advection, January 2016', fontsize=28)
+    # Salinity
+    ax = subplot(gs2[0,1])
+    pcolor(lat, z, upwind_salt, vmin=salt_min, vmax=salt_max, cmap='jet')
+    contour(lat, z, upwind_salt, levels=[salt_contour], color='black')
+    xlim([lat_min, lat_max])
+    ylim([depth_min, depth_max])
+    ax.set_xticks(lat_ticks)
+    ax.set_xticklabels(lat_labels, fontsize=16)
+    ax.set_yticks(depth_ticks)
+    ax.set_yticklabels([])
+    # Akima advection
+    gs3 = GridSpec(1,2)
+    gs3.update(left=0.1, right=0.95, bottom=0.2925, top=0.475, wspace=0.08)
+    # Temperature
+    ax = subplot(gs3[0,0])
+    pcolor(lat, z, akima_temp, vmin=temp_min, vmax=temp_max, cmap='jet')
+    contour(lat, z, akima_temp, levels=[temp_contour], color='black')
+    ylabel('Depth (m)', fontsize=18)
+    xlim([lat_min, lat_max])
+    ylim([depth_min, depth_max])
+    ax.set_xticks(lat_ticks)
+    ax.set_xticklabels(lat_labels, fontsize=16)
+    ax.set_yticks(depth_ticks)
+    ax.set_yticklabels(depth_labels, fontsize=16)
+    text(-52, 300, 'c) Akima advection, January 2016', fontsize=28)
+    # Salinity
+    ax = subplot(gs3[0,1])
+    pcolor(lat, z, akima_salt, vmin=salt_min, vmax=salt_max, cmap='jet')
+    contour(lat, z, akima_salt, levels=[salt_contour], color='black')
+    xlim([lat_min, lat_max])
+    ylim([depth_min, depth_max])
+    ax.set_xticks(lat_ticks)
+    ax.set_xticklabels(lat_labels, fontsize=16)
+    ax.set_yticks(depth_ticks)
+    ax.set_yticklabels([])
+    # Split advection
+    gs4 = GridSpec(1,2)
+    gs4.update(left=0.1, right=0.95, bottom=0.06, top=0.2425, wspace=0.08)
+    # Temperature
+    ax = subplot(gs4[0,0])
+    img = pcolor(lat, z, split_temp, vmin=temp_min, vmax=temp_max, cmap='jet')
+    contour(lat, z, split_temp, levels=[temp_contour], color='black')
+    ylabel('Depth (m)', fontsize=18)
+    xlim([lat_min, lat_max])
+    ylim([depth_min, depth_max])
+    ax.set_xticks(lat_ticks)
+    ax.set_xticklabels(lat_labels, fontsize=16)
+    ax.set_yticks(depth_ticks)
+    ax.set_yticklabels(depth_labels, fontsize=16)
+    text(-52, 300, 'd) RSUP3 advection, January 2016', fontsize=28)
+    # Add a colorbar for temperature
+    cbaxes = fig.add_axes([0.17, 0.015, 0.3, 0.015])
+    cbar = colorbar(img, orientation='horizontal', cax=cbaxes, extend='both', ticks=arange(temp_min, temp_max+2, 2))
+    cbar.ax.tick_params(labelsize=16)
+    # Salinity
+    ax = subplot(gs4[0,1])
+    img = pcolor(lat, z, split_salt, vmin=salt_min, vmax=salt_max, cmap='jet')
+    contour(lat, z, split_salt, levels=[salt_contour], color='black')
+    xlim([lat_min, lat_max])
+    ylim([depth_min, depth_max])
+    ax.set_xticks(lat_ticks)
+    ax.set_xticklabels(lat_labels, fontsize=16)
+    ax.set_yticks(depth_ticks)
+    ax.set_yticklabels([])
+    # Add a colorbar for salinity
+    cbaxes = fig.add_axes([0.6, 0.015, 0.3, 0.02])
+    cbar = colorbar(img, orientation='horizontal', cax=cbaxes, extend='both', ticks=arange(salt_min+0.1, salt_max+0.1, 0.2))
+    cbar.ax.tick_params(labelsize=16)
+    fig.show()
+    fig.savefig('bugs_drift.png')
+
+
+# Command-line interface
+if __name__ == "__main__":
+
+    grid_file = raw_input("Path to ROMS grid file: ")
+    upwind_file = raw_input("Path to upwind advection file containing monthly averaged temperature and salinity for January 2016: ")
+    akima_file = raw_input("Path to Akima advection file containing monthly averaged temperature and salinity for January 2016: ")
+    split_file = raw_input("Path to split advection file containing monthly averaged temperature and salinity for January 2016: ")
+    bugs_drift_slices(grid_file, upwind_file, akima_file, split_file)
diff --git a/mip_aice_minmax_nsidc.py b/mip_aice_minmax_nsidc.py
index c36e8e2..301fec2 100644
--- a/mip_aice_minmax_nsidc.py
+++ b/mip_aice_minmax_nsidc.py
@@ -17,12 +17,15 @@
 # cice_file = path to CICE output file containing 5-day averages for the entire
 #             simulation
 # cice_log = path to CICE logfile from timeseries_seaice_extent.py
-# fesom_mesh_path = path to FESOM mesh directory
-# fesom_output_dir = path to FESOM output directory containing one ice.mean.nc
-#                    file for each year (5-day averages)
-# fesom_log = path to FESOM logfile from timeseries_seaice_extent.py in
-#             fesomtools
-def mip_aice_minmax_nsidc (cice_file, cice_log, fesom_mesh_path, fesom_output_dir, fesom_log):
+# fesom_mesh_path_lr, fesom_mesh_path_hr = paths to FESOM mesh directories for
+#                     low-res and high-res mesh respectively
+# fesom_output_dir_lr, fesom_output_dir_hr = paths to FESOM output directories
+#                      containing one ice.mean.nc file for each year (5-day
+#                      averages), low-res and high-res respectively
+# fesom_log_lr, fesom_log_hr = paths to FESOM logfiles from
+#               timeseries_seaice_extent.py in fesomtools, low-res and high-res
+#               respectively
+def mip_aice_minmax_nsidc (cice_file, cice_log, fesom_mesh_path_lr, fesom_output_dir_lr, fesom_log_lr, fesom_mesh_path_hr, fesom_output_dir_hr, fesom_log_hr):
 
     # Range of years to process (exclude 2016 because no NSIDC data)
     start_year = 1992
@@ -247,60 +250,94 @@ def mip_aice_minmax_nsidc (cice_file, cice_log, fesom_mesh_path, fesom_output_di
         cice_feb_extent.append(cice_extent[(year-start_year)*12+1])
         cice_sep_extent.append(cice_extent[(year-start_year)*12+8])
 
-    print 'Processing FESOM'
+    print 'Processing FESOM low-res'
 
     # First build the grid
-    elements, patches = make_patches(fesom_mesh_path, circumpolar, mask_cavities)
+    elements_lr, patches_lr = make_patches(fesom_mesh_path_lr, circumpolar, mask_cavities)
     # Get averages for February and September
     # Start with first year just to initialise the arrays with the right size
     print '...monthly average for ' + str(start_year)
-    fesom_feb_nodes = monthly_avg(fesom_output_dir + fesom_file_head + '.' + str(start_year) + '.ice.mean.nc', 'area', 1)
-    fesom_sep_nodes = monthly_avg(fesom_output_dir + fesom_file_head + '.' + str(start_year) + '.ice.mean.nc', 'area', 8)
+    fesom_feb_nodes_lr = monthly_avg(fesom_output_dir_lr + fesom_file_head + '.' + str(start_year) + '.ice.mean.nc', 'area', 1)
+    fesom_sep_nodes_lr = monthly_avg(fesom_output_dir_lr + fesom_file_head + '.' + str(start_year) + '.ice.mean.nc', 'area', 8)
     # Loop over the rest of the years
     for year in range(start_year+1, end_year+1):
         print '...monthly average for ' + str(year)
-        fesom_feb_nodes = fesom_feb_nodes + monthly_avg(fesom_output_dir + fesom_file_head + '.' + str(year) + '.ice.mean.nc', 'area', 1)
-        fesom_sep_nodes = fesom_sep_nodes + monthly_avg(fesom_output_dir + fesom_file_head + '.' + str(year) + '.ice.mean.nc', 'area', 8)
+        fesom_feb_nodes_lr = fesom_feb_nodes_lr + monthly_avg(fesom_output_dir_lr + fesom_file_head + '.' + str(year) + '.ice.mean.nc', 'area', 1)
+        fesom_sep_nodes_lr = fesom_sep_nodes_lr + monthly_avg(fesom_output_dir_lr + fesom_file_head + '.' + str(year) + '.ice.mean.nc', 'area', 8)
     # Convert from integrals to averages
-    fesom_feb_nodes = fesom_feb_nodes/num_years
-    fesom_sep_nodes = fesom_sep_nodes/num_years
+    fesom_feb_nodes_lr = fesom_feb_nodes_lr/num_years
+    fesom_sep_nodes_lr = fesom_sep_nodes_lr/num_years
     # Find element-averages
-    fesom_feb = []
-    fesom_sep = []
-    for elm in elements:
+    fesom_feb_lr = []
+    fesom_sep_lr = []
+    for elm in elements_lr:
         if not elm.cavity:
             # Average over 3 component nodes
-            fesom_feb.append(mean(array([fesom_feb_nodes[elm.nodes[0].id], fesom_feb_nodes[elm.nodes[1].id], fesom_feb_nodes[elm.nodes[2].id]])))
-            fesom_sep.append(mean(array([fesom_sep_nodes[elm.nodes[0].id], fesom_sep_nodes[elm.nodes[1].id], fesom_sep_nodes[elm.nodes[2].id]])))
-    fesom_feb = array(fesom_feb)
-    fesom_sep = array(fesom_sep)
+            fesom_feb_lr.append(mean(array([fesom_feb_nodes_lr[elm.nodes[0].id], fesom_feb_nodes_lr[elm.nodes[1].id], fesom_feb_nodes_lr[elm.nodes[2].id]])))
+            fesom_sep_lr.append(mean(array([fesom_sep_nodes_lr[elm.nodes[0].id], fesom_sep_nodes_lr[elm.nodes[1].id], fesom_sep_nodes_lr[elm.nodes[2].id]])))
+    fesom_feb_lr = array(fesom_feb_lr)
+    fesom_sep_lr = array(fesom_sep_lr)
 
     # Get extent timeseries
     # Read 5-day logfile
-    fesom_extent_5day = []
-    f = open(fesom_log, 'r')    
+    fesom_extent_5day_lr = []
+    f = open(fesom_log_lr, 'r')    
     f.readline()
     for line in f:
-        fesom_extent_5day.append(float(line))
+        fesom_extent_5day_lr.append(float(line))
     f.close()
     # Initialise monthly arrays
-    fesom_feb_extent = zeros(num_years)
-    fesom_sep_extent = zeros(num_years)
+    fesom_feb_extent_lr = zeros(num_years)
+    fesom_sep_extent_lr = zeros(num_years)
     for year in range(start_year, end_year+1):
         # First timestep of year in 5-day logfile
         t0 = (year-start_year)*73
         # Monthly averages are hard-coded and ugly
         # Feburary: 4/5 of index 7, indices 8-11, and 4/5 of index 12
-        fesom_feb_extent[year-start_year] = (fesom_extent_5day[t0+6]*4 + sum(fesom_extent_5day[t0+7:t0+11]*5) + fesom_extent_5day[t0+11]*4)/28.0
+        fesom_feb_extent_lr[year-start_year] = (fesom_extent_5day_lr[t0+6]*4 + sum(fesom_extent_5day_lr[t0+7:t0+11]*5) + fesom_extent_5day_lr[t0+11]*4)/28.0
         # September: 2/5 of index 49, indices 50-54, 3/5 of index 55
-        fesom_sep_extent[year-start_year] = (fesom_extent_5day[t0+48]*2 + sum(fesom_extent_5day[t0+49:t0+54]*5) + fesom_extent_5day[t0+54]*3)/30.0
+        fesom_sep_extent_lr[year-start_year] = (fesom_extent_5day_lr[t0+48]*2 + sum(fesom_extent_5day_lr[t0+49:t0+54]*5) + fesom_extent_5day_lr[t0+54]*3)/30.0
+
+    print 'Processing FESOM high-res'
+
+    elements_hr, patches_hr = make_patches(fesom_mesh_path_hr, circumpolar, mask_cavities)
+    print '...monthly average for ' + str(start_year)
+    fesom_feb_nodes_hr = monthly_avg(fesom_output_dir_hr + fesom_file_head + '.' + str(start_year) + '.ice.mean.nc', 'area', 1)
+    fesom_sep_nodes_hr = monthly_avg(fesom_output_dir_hr + fesom_file_head + '.' + str(start_year) + '.ice.mean.nc', 'area', 8)
+    for year in range(start_year+1, end_year+1):
+        print '...monthly average for ' + str(year)
+        fesom_feb_nodes_hr = fesom_feb_nodes_hr + monthly_avg(fesom_output_dir_hr + fesom_file_head + '.' + str(year) + '.ice.mean.nc', 'area', 1)
+        fesom_sep_nodes_hr = fesom_sep_nodes_hr + monthly_avg(fesom_output_dir_hr + fesom_file_head + '.' + str(year) + '.ice.mean.nc', 'area', 8)
+    fesom_feb_nodes_hr = fesom_feb_nodes_hr/num_years
+    fesom_sep_nodes_hr = fesom_sep_nodes_hr/num_years
+    fesom_feb_hr = []
+    fesom_sep_hr = []
+    for elm in elements_hr:
+        if not elm.cavity:
+            fesom_feb_hr.append(mean(array([fesom_feb_nodes_hr[elm.nodes[0].id], fesom_feb_nodes_hr[elm.nodes[1].id], fesom_feb_nodes_hr[elm.nodes[2].id]])))
+            fesom_sep_hr.append(mean(array([fesom_sep_nodes_hr[elm.nodes[0].id], fesom_sep_nodes_hr[elm.nodes[1].id], fesom_sep_nodes_hr[elm.nodes[2].id]])))
+    fesom_feb_hr = array(fesom_feb_hr)
+    fesom_sep_hr = array(fesom_sep_hr)
+
+    fesom_extent_5day_hr = []
+    f = open(fesom_log_hr, 'r')    
+    f.readline()
+    for line in f:
+        fesom_extent_5day_hr.append(float(line))
+    f.close()
+    fesom_feb_extent_hr = zeros(num_years)
+    fesom_sep_extent_hr = zeros(num_years)
+    for year in range(start_year, end_year+1):
+        t0 = (year-start_year)*73
+        fesom_feb_extent_hr[year-start_year] = (fesom_extent_5day_hr[t0+6]*4 + sum(fesom_extent_5day_hr[t0+7:t0+11]*5) + fesom_extent_5day_hr[t0+11]*4)/28.0
+        fesom_sep_extent_hr[year-start_year] = (fesom_extent_5day_hr[t0+48]*2 + sum(fesom_extent_5day_hr[t0+49:t0+54]*5) + fesom_extent_5day_hr[t0+54]*3)/30.0
 
     time_axis = arange(start_year, end_year+1)
 
     print 'Plotting'
-    fig = figure(figsize=(20,10))
-    gs1 = GridSpec(2, 3)
-    gs1.update(left=0.12, right=0.7, wspace=0.05, hspace=0.05)
+    fig = figure(figsize=(24,10))
+    gs1 = GridSpec(2, 4)
+    gs1.update(left=0.1, right=0.77, wspace=0.04, hspace=0.05)
     # NSIDC, February
     ax = subplot(gs1[0, 0], aspect='equal')
     img = pcolor(nsidc_x, nsidc_y, nsidc_feb, vmin=0, vmax=1, cmap='jet')
@@ -318,10 +355,22 @@ def mip_aice_minmax_nsidc (cice_file, cice_log, fesom_mesh_path, fesom_output_di
     ax.set_xticks([])
     ax.set_yticks([])
     title('MetROMS', fontsize=24)
-    # FESOM, February
+    # FESOM low-res, February
     ax = subplot(gs1[0, 2], aspect='equal')
-    img = PatchCollection(patches, cmap='jet')
-    img.set_array(fesom_feb)
+    img = PatchCollection(patches_lr, cmap='jet')
+    img.set_array(fesom_feb_lr)
+    img.set_clim(vmin=0, vmax=1)
+    img.set_edgecolor('face')
+    ax.add_collection(img)
+    xlim([bdry1, bdry2])
+    ylim([bdry3, bdry4])
+    ax.set_xticks([])
+    ax.set_yticks([])
+    title('FESOM (low-res)', fontsize=24)
+    # FESOM high-res, February
+    ax = subplot(gs1[0, 3], aspect='equal')
+    img = PatchCollection(patches_hr, cmap='jet')
+    img.set_array(fesom_feb_hr)
     img.set_clim(vmin=0, vmax=1)
     img.set_edgecolor('face')
     ax.add_collection(img)
@@ -347,10 +396,21 @@ def mip_aice_minmax_nsidc (cice_file, cice_log, fesom_mesh_path, fesom_output_di
     ylim([bdry3, bdry4])
     ax.set_xticks([])
     ax.set_yticks([])    
-    # FESOM, September
+    # FESOM low-res, September
     ax = subplot(gs1[1, 2], aspect='equal')
-    img = PatchCollection(patches, cmap='jet')
-    img.set_array(fesom_sep)
+    img = PatchCollection(patches_lr, cmap='jet')
+    img.set_array(fesom_sep_lr)
+    img.set_clim(vmin=0, vmax=1)
+    img.set_edgecolor('face')
+    ax.add_collection(img)
+    xlim([bdry1, bdry2])
+    ylim([bdry3, bdry4])
+    ax.set_xticks([])
+    ax.set_yticks([])
+    # FESOM high-res, September
+    ax = subplot(gs1[1, 3], aspect='equal')
+    img = PatchCollection(patches_hr, cmap='jet')
+    img.set_array(fesom_sep_hr)
     img.set_clim(vmin=0, vmax=1)
     img.set_edgecolor('face')
     ax.add_collection(img)
@@ -359,17 +419,18 @@ def mip_aice_minmax_nsidc (cice_file, cice_log, fesom_mesh_path, fesom_output_di
     ax.set_xticks([])
     ax.set_yticks([])
     # Add a colourbar at the bottom
-    cbaxes = fig.add_axes([0.12, 0.04, 0.3, 0.04])
-    cbar = colorbar(img, orientation='horizontal', ticks=arange(0,1+0.25,0.25), cax=cbaxes, extend='min')
+    cbaxes = fig.add_axes([0.1, 0.04, 0.25, 0.04])
+    cbar = colorbar(img, orientation='horizontal', ticks=arange(0,1+0.25,0.25), cax=cbaxes)
     cbar.ax.tick_params(labelsize=20)
     # Add extent timeseries on rightmost column, with more space for labels
     gs2 = GridSpec(2, 1)
-    gs2.update(left=0.73, right=0.95, wspace=0.1, hspace=0.15)
+    gs2.update(left=0.79, right=0.95, hspace=0.15)
     # February
     ax = subplot(gs2[0, 0])
     ax.plot(time_axis, nsidc_feb_extent, color='black', linewidth=2, linestyle='dashed')
     ax.plot(time_axis, cice_feb_extent, color='blue', linewidth=1.5)
-    ax.plot(time_axis, fesom_feb_extent, color='green', linewidth=1.5)
+    ax.plot(time_axis, fesom_feb_extent_lr, color='green', linewidth=1.5)
+    ax.plot(time_axis, fesom_feb_extent_hr, color='magenta', linewidth=1.5)
     xlim([start_year, end_year])
     ax.set_yticks(arange(0,4+0.5,0.5))
     ax.set_yticklabels(['0', '', '1', '', '2', '', '3', '', '4'])
@@ -381,14 +442,15 @@ def mip_aice_minmax_nsidc (cice_file, cice_log, fesom_mesh_path, fesom_output_di
     ax = subplot(gs2[1, 0])
     ax.plot(time_axis, nsidc_sep_extent, color='black', label='NSIDC', linewidth=2, linestyle='dashed')
     ax.plot(time_axis, cice_sep_extent, color='blue', label='MetROMS', linewidth=1.5)
-    ax.plot(time_axis, fesom_sep_extent, color='green', label='FESOM (high-res)', linewidth=1.5)
+    ax.plot(time_axis, fesom_sep_extent_lr, color='green', label='FESOM (low-res)', linewidth=1.5)
+    ax.plot(time_axis, fesom_sep_extent_hr, color='magenta', label='FESOM (high-res)', linewidth=1.5)
     xlim([start_year, end_year])
     ax.set_yticks(arange(16,23+1,1))
     ax.set_yticklabels(['16', '', '18', '', '20', '', '22', ''])
     ax.tick_params(axis='x', labelsize=20)
     ax.tick_params(axis='y', labelsize=20)
     grid(True)
-    ax.legend(bbox_to_anchor=(1.04,-0.07), ncol=3, fontsize=20)
+    ax.legend(bbox_to_anchor=(1.04,-0.07), ncol=4, fontsize=20)
     fig.show()
     fig.savefig('aice_minmax_nsidc.png')
 
@@ -398,10 +460,13 @@ def mip_aice_minmax_nsidc (cice_file, cice_log, fesom_mesh_path, fesom_output_di
 
     cice_file = raw_input("Path to CICE output file containing data for the entire simulation: ")
     cice_log = raw_input("Path to CICE sea ice extent logfile: ")
-    fesom_mesh_path = raw_input("Path to FESOM mesh directory: ")
-    fesom_output_dir = raw_input("Path to FESOM output directory containing one ice.mean.nc file for each year: ")
-    fesom_log = raw_input("Path to FESOM sea ice extent logfile: ")
-    mip_aice_minmax_nsidc(cice_file, cice_log, fesom_mesh_path, fesom_output_dir, fesom_log)
+    fesom_mesh_path_lr = raw_input("Path to FESOM low-res mesh directory: ")
+    fesom_output_dir_lr = raw_input("Path to FESOM low-res output directory containing one ice.mean.nc file for each year: ")
+    fesom_log_lr = raw_input("Path to FESOM low-res sea ice extent logfile: ")
+    fesom_mesh_path_hr = raw_input("Path to FESOM high-res mesh directory: ")
+    fesom_output_dir_hr = raw_input("Path to FESOM high-res output directory containing one ice.mean.nc file for each year: ")
+    fesom_log_hr = raw_input("Path to FESOM high-res sea ice extent logfile: ")
+    mip_aice_minmax_nsidc(cice_file, cice_log, fesom_mesh_path_lr, fesom_output_dir_lr, fesom_log_lr, fesom_mesh_path_hr, fesom_output_dir_hr, fesom_log_hr)
     
          
     
diff --git a/mip_calc_massloss.py b/mip_calc_massloss.py
index f757fe3..cc9b0fa 100644
--- a/mip_calc_massloss.py
+++ b/mip_calc_massloss.py
@@ -5,24 +5,34 @@
 # estimates for 2003-2008.
 # Input:
 # roms_logfile = path to ROMS logfile from timeseries_massloss.py
+# roms_logfile_bs = path to ROMS logfile from
+#                   timeseries_massloss_bellingshausen.py
 # fesom_logfile_lr, fesom_logfile_hr = paths to FESOM logfiles from
 #                   timeseries_massloss.py in the fesomtools repository, for
 #                   low-res and high-res respectively
-def mip_calc_massloss (roms_logfile, fesom_logfile_lr, fesom_logfile_hr):
+# fesom_logfile_bs_lr, fesom_logfile_bs_hr = paths to FESOM logfiles from
+#                      timeseries_massloss_bellingshausen.py in the fesomtools
+#                      repository, for low-res and high-res respectively
+def mip_calc_massloss (roms_logfile, roms_logfile_bs, fesom_logfile_lr, fesom_logfile_bs_lr, fesom_logfile_hr, fesom_logfile_bs_hr):
 
     # Year simulations start
     year_start = 1992
-    # Years to avearge over
+    # Years to average over
     calc_start = 2002
     calc_end = 2016
     # Number of output steps per year in FESOM
     peryear = 365/5
     # Name of each ice shelf
-    names = ['Total Mass Loss', 'Larsen D Ice Shelf', 'Larsen C Ice Shelf', 'Wilkins & George VI & Stange Ice Shelves', 'Ronne-Filchner Ice Shelf', 'Abbot Ice Shelf', 'Pine Island Glacier Ice Shelf', 'Thwaites Ice Shelf', 'Dotson Ice Shelf', 'Getz Ice Shelf', 'Nickerson Ice Shelf', 'Sulzberger Ice Shelf', 'Mertz Ice Shelf', 'Totten & Moscow University Ice Shelves', 'Shackleton Ice Shelf', 'West Ice Shelf', 'Amery Ice Shelf', 'Prince Harald Ice Shelf', 'Baudouin & Borchgrevink Ice Shelves', 'Lazarev Ice Shelf', 'Nivl Ice Shelf', 'Fimbul & Jelbart & Ekstrom Ice Shelves', 'Brunt & Riiser-Larsen Ice Shelves', 'Ross Ice Shelf']
+    names = ['Total Mass Loss', 'Larsen D Ice Shelf', 'Larsen C Ice Shelf', 'Wilkins & George VI & Stange & Bach Ice Shelves', 'Ronne-Filchner Ice Shelf', 'Abbot Ice Shelf', 'Pine Island Glacier Ice Shelf', 'Thwaites Ice Shelf', 'Dotson Ice Shelf', 'Getz Ice Shelf', 'Nickerson Ice Shelf', 'Sulzberger Ice Shelf', 'Mertz Ice Shelf', 'Totten & Moscow University Ice Shelves', 'Shackleton Ice Shelf', 'West Ice Shelf', 'Amery Ice Shelf', 'Prince Harald Ice Shelf', 'Baudouin & Borchgrevink Ice Shelves', 'Lazarev Ice Shelf', 'Nivl Ice Shelf', 'Fimbul & Jelbart & Ekstrom Ice Shelves', 'Brunt & Riiser-Larsen Ice Shelves', 'Ross Ice Shelf']
     # Observed mass loss (Rignot 2013) and uncertainty for each ice shelf, in Gt/y
     obs_massloss = [1325, 1.4, 20.7, 135.4, 155.4, 51.8, 101.2, 97.5, 45.2, 144.9, 4.2, 18.2, 7.9, 90.6, 72.6, 27.2, 35.5, -2, 21.6, 6.3, 3.9, 26.8, 9.7, 47.7]
     obs_massloss_error = [235, 14, 67, 40, 45, 19, 8, 7, 4, 14, 2, 3, 3, 8, 15, 10, 23, 3, 18, 2, 2, 14, 16, 34]
     num_shelves = len(obs_massloss)
+    # Some Bellingshausen ice shelves were split up later
+    names_bs = ['Wilkins Ice Shelf', 'Stange Ice Shelf', 'George VI Ice Shelf']
+    obs_massloss_bs = [18.4, 28, 89]
+    obs_massloss_error_bs = [17, 6, 17]
+    num_shelves_bs = len(obs_massloss_bs)
 
     # Read ROMS logfile
     roms_time = []
@@ -60,6 +70,35 @@ def mip_calc_massloss (roms_logfile, fesom_logfile_lr, fesom_logfile_hr):
         t_end = nonzero(roms_time >= calc_end+1)[0][0]
     roms_massloss = mean(roms_massloss_ts[:,t_start:t_end], axis=1)
 
+    # Repeat for Bellingshausen
+    f = open(roms_logfile_bs, 'r')
+    f.readline()
+    # Skip the time values (should be the same)
+    for line in f:
+        try:
+            tmp = float(line)
+        except(ValueError):
+            # Reached the header for the next variable
+            break       
+    roms_massloss_bs_ts = empty([num_shelves_bs, len(roms_time)])
+    index = 0
+    while index < num_shelves_bs:
+        t = 0
+        for line in f:
+            try:
+                roms_massloss_bs_ts[index, t] = float(line)
+                t += 1
+            except(ValueError):
+                break
+        index +=1
+    f.close()
+    t_start = nonzero(roms_time >= calc_start)[0][0]
+    if calc_end == 2016:
+        t_end = size(roms_time)
+    else:
+        t_end = nonzero(roms_time >= calc_end+1)[0][0]
+    roms_massloss_bs = mean(roms_massloss_bs_ts[:,t_start:t_end], axis=1)
+
     # Read FESOM timeseries
     # Low-res
     tmp = []
@@ -94,6 +133,24 @@ def mip_calc_massloss (roms_logfile, fesom_logfile_lr, fesom_logfile_hr):
     f.close()
     # Average between given years
     fesom_massloss_lr = mean(fesom_massloss_ts_lr[:,peryear*(calc_start-year_start):peryear*(calc_end+1-year_start)], axis=1)
+
+    # Repeat for Bellingshausen
+    f = open(fesom_logfile_bs_lr, 'r')
+    f.readline()
+    fesom_massloss_bs_ts_lr = empty([num_shelves_bs, num_time])
+    index = 0
+    while index < num_shelves_bs:
+        t = 0
+        for line in f:
+            try:
+                fesom_massloss_bs_ts_lr[index,t] = float(line)
+                t += 1
+            except(ValueError):
+                break
+        index += 1
+    f.close()
+    fesom_massloss_bs_lr = mean(fesom_massloss_bs_ts_lr[:,peryear*(calc_start-year_start):peryear*(calc_end+1-year_start)], axis=1)
+
     # Repeat for high-res
     tmp = []
     f = open(fesom_logfile_hr, 'r')
@@ -120,6 +177,23 @@ def mip_calc_massloss (roms_logfile, fesom_logfile_lr, fesom_logfile_hr):
     f.close()
     fesom_massloss_hr = mean(fesom_massloss_ts_hr[:,peryear*(calc_start-year_start):peryear*(calc_end+1-year_start)], axis=1)
 
+    # High-res Bellingshausen
+    f = open(fesom_logfile_bs_hr, 'r')
+    f.readline()
+    fesom_massloss_bs_ts_hr = empty([num_shelves_bs, num_time])
+    index = 0
+    while index < num_shelves_bs:
+        t = 0
+        for line in f:
+            try:
+                fesom_massloss_bs_ts_hr[index,t] = float(line)
+                t += 1
+            except(ValueError):
+                break
+        index += 1
+    f.close()
+    fesom_massloss_bs_hr = mean(fesom_massloss_bs_ts_hr[:,peryear*(calc_start-year_start):peryear*(calc_end+1-year_start)], axis=1)
+
     # Loop over ice shelves
     for index in range(num_shelves):
         print names[index]
@@ -131,14 +205,28 @@ def mip_calc_massloss (roms_logfile, fesom_logfile_lr, fesom_logfile_hr):
         massloss_high = obs_massloss[index] + obs_massloss_error[index]
         print 'Rignot: ' + str(massloss_low) + '-' + str(massloss_high)
 
+    # Repeat for Bellingshausen ice shelves
+    for index in range(num_shelves_bs):
+        print names_bs[index]
+        print 'MetROMS: ' + str(roms_massloss_bs[index])
+        print 'FESOM low-res: ' + str(fesom_massloss_bs_lr[index])
+        print 'FESOM high-res: ' + str(fesom_massloss_bs_hr[index])
+        # Find the range of observations
+        massloss_low = obs_massloss_bs[index] - obs_massloss_error_bs[index]
+        massloss_high = obs_massloss_bs[index] + obs_massloss_error_bs[index]
+        print 'Rignot: ' + str(massloss_low) + '-' + str(massloss_high)
+
 
 # Command-line interface
 if __name__ == "__main__":
 
     roms_logfile = raw_input("Path to ROMS logfile from timeseries_massloss.py: ")
+    roms_logfile_bs = raw_input("Path to ROMS logfile from timeseries_massloss_bellingshausen.py: ")
     fesom_logfile_lr = raw_input("Path to FESOM low-res logfile from timeseries_massloss.py: ")
+    fesom_logfile_bs_lr = raw_input("Path to FESOM low-res logfile from timeseries_massloss_bellingshausen.py: ")
     fesom_logfile_hr = raw_input("Path to FESOM high-res logfile from timeseries_massloss.py: ")
-    mip_calc_massloss(roms_logfile, fesom_logfile_lr, fesom_logfile_hr)
+    fesom_logfile_bs_hr = raw_input("Path to FESOM high-res logfile from timeseries_massloss_bellingshausen.py: ")
+    mip_calc_massloss(roms_logfile, roms_logfile_bs, fesom_logfile_lr, fesom_logfile_bs_lr, fesom_logfile_hr, fesom_logfile_bs_hr)
 
 
 
diff --git a/mip_drift_slices.py b/mip_drift_slices.py
index 50f53a1..ac09f4f 100644
--- a/mip_drift_slices.py
+++ b/mip_drift_slices.py
@@ -22,10 +22,12 @@
 # roms_grid = path to ROMS grid file
 # roms_file = path to file containing Jan 2016 monthly average of temperature
 #             and salinity in ROMS
-# fesom_mesh_path = path to FESOM mesh directory
-# fesom_file = path to file containing Jan 2016 monthly average of temperature
-#              and salinity in FESOM
-def mip_drift_slices (roms_grid, roms_file, fesom_mesh_path, fesom_file):
+# fesom_mesh_path_lr, fesom_mesh_path_hr = paths to FESOM mesh directories for
+#                     low-res and high-res respectively
+# fesom_file_lr, fesom_file_hr = paths to files containing Jan 2016 monthly
+#                averages of temperature and salinity, in low-res FESOM and
+#                high-res FESOM respectively
+def mip_drift_slices (roms_grid, roms_file, fesom_mesh_path_lr, fesom_file_lr, fesom_mesh_path_hr, fesom_file_hr):
 
     # Paths to ECCO2 files with initial conditions for temp and salt
     ecco_temp_file = '/short/m68/kaa561/metroms_iceshelf/data/originals/ECCO2/THETA.1440x720x50.199201.nc'
@@ -51,8 +53,8 @@ def mip_drift_slices (roms_grid, roms_file, fesom_mesh_path, fesom_file):
     temp_contour = 0.75
     salt_contour = 34.5
     # Parameters for FESOM regular grid interpolation (needed for contours)
-    num_lat = 1000
-    num_depth = 500
+    num_lat = 500
+    num_depth = 250
     r = 6.371e6
     deg2rad = pi/180.0
 
@@ -104,42 +106,42 @@ def mip_drift_slices (roms_grid, roms_file, fesom_mesh_path, fesom_file):
     if lon0 > 180:
         lon0 -= 360
 
-    print 'Processing FESOM'
+    print 'Processing low-res FESOM'
     # Build regular elements
-    elements = fesom_grid(fesom_mesh_path)
+    elements_lr = fesom_grid(fesom_mesh_path_lr)
     # Read temperature and salinity
-    id = Dataset(fesom_file, 'r')
-    fesom_temp_nodes = id.variables['temp'][0,:]
-    fesom_salt_nodes = id.variables['salt'][0,:]
+    id = Dataset(fesom_file_lr, 'r')
+    fesom_temp_nodes_lr = id.variables['temp'][0,:]
+    fesom_salt_nodes_lr = id.variables['salt'][0,:]
     id.close()
     # Make SideElements
-    selements_temp = fesom_sidegrid(elements, fesom_temp_nodes, lon0, lat_max)
-    selements_salt = fesom_sidegrid(elements, fesom_salt_nodes, lon0, lat_max)
+    selements_temp_lr = fesom_sidegrid(elements_lr, fesom_temp_nodes_lr, lon0, lat_max)
+    selements_salt_lr = fesom_sidegrid(elements_lr, fesom_salt_nodes_lr, lon0, lat_max)
     # Build an array of quadrilateral patches for the plot, and of data values
     # corresponding to each SideElement
-    patches = []
-    fesom_temp = []
-    for selm in selements_temp:
+    patches_lr = []
+    fesom_temp_lr = []
+    for selm in selements_temp_lr:
         # Make patch
         coord = transpose(vstack((selm.y, selm.z)))
-        patches.append(Polygon(coord, True, linewidth=0.))
+        patches_lr.append(Polygon(coord, True, linewidth=0.))
         # Save data value
-        fesom_temp.append(selm.var)
-    # Repeat for the other variables
-    fesom_salt = []
-    for selm in selements_salt:
-        fesom_salt.append(selm.var)
+        fesom_temp_lr.append(selm.var)
+    # Repeat for salinity
+    fesom_salt_lr = []
+    for selm in selements_salt_lr:
+        fesom_salt_lr.append(selm.var)
     # Interpolate to regular grid so we can overlay contours
     lat_reg = linspace(lat_min, lat_max, num_lat)
     depth_reg = linspace(-depth_max, -depth_min, num_depth)
-    temp_reg = zeros([num_depth, num_lat])
-    salt_reg = zeros([num_depth, num_lat])
-    temp_reg[:,:] = NaN
-    salt_reg[:,:] = NaN
+    temp_reg_lr = zeros([num_depth, num_lat])
+    salt_reg_lr = zeros([num_depth, num_lat])
+    temp_reg_lr[:,:] = NaN
+    salt_reg_lr[:,:] = NaN
     # For each element, check if a point on the regular grid lies
     # within. If so, do barycentric interpolation to that point, at each
     # depth on the regular grid.
-    for elm in elements:
+    for elm in elements_lr:
         # Check if this element crosses lon0
         if amin(elm.lon) < lon0 and amax(elm.lon) > lon0:
             # Check if we are within the latitude bounds
@@ -185,17 +187,80 @@ def mip_drift_slices (roms_grid, roms_file, fesom_mesh_path, fesom_file):
                                     node_vals_temp.append(NaN)
                                     node_vals_salt.append(NaN)
                                 else:
-                                    node_vals_temp.append(coeff1*fesom_temp_nodes[id1] + coeff2*fesom_temp_nodes[id2])
-                                    node_vals_salt.append(coeff1*fesom_salt_nodes[id1] + coeff2*fesom_salt_nodes[id2])
+                                    node_vals_temp.append(coeff1*fesom_temp_nodes_lr[id1] + coeff2*fesom_temp_nodes_lr[id2])
+                                    node_vals_salt.append(coeff1*fesom_salt_nodes_lr[id1] + coeff2*fesom_salt_nodes_lr[id2])
                             if any(isnan(node_vals_temp)):
                                 pass
                             else:
                                 # Barycentric interpolation for the value at
                                 # lon0, lat0
-                                temp_reg[k,j] = sum(array(cff)*array(node_vals_temp))
-                                salt_reg[k,j] = sum(array(cff)*array(node_vals_salt))
-    temp_reg = ma.masked_where(isnan(temp_reg), temp_reg)
-    salt_reg = ma.masked_where(isnan(salt_reg), salt_reg)
+                                temp_reg_lr[k,j] = sum(array(cff)*array(node_vals_temp))
+                                salt_reg_lr[k,j] = sum(array(cff)*array(node_vals_salt))
+    temp_reg_lr = ma.masked_where(isnan(temp_reg_lr), temp_reg_lr)
+    salt_reg_lr = ma.masked_where(isnan(salt_reg_lr), salt_reg_lr)
+
+    print 'Processing high-res FESOM'
+    elements_hr = fesom_grid(fesom_mesh_path_hr)
+    id = Dataset(fesom_file_hr, 'r')
+    fesom_temp_nodes_hr = id.variables['temp'][0,:]
+    fesom_salt_nodes_hr = id.variables['salt'][0,:]
+    id.close()
+    selements_temp_hr = fesom_sidegrid(elements_hr, fesom_temp_nodes_hr, lon0, lat_max)
+    selements_salt_hr = fesom_sidegrid(elements_hr, fesom_salt_nodes_hr, lon0, lat_max)
+    patches_hr = []
+    fesom_temp_hr = []
+    for selm in selements_temp_hr:
+        coord = transpose(vstack((selm.y, selm.z)))
+        patches_hr.append(Polygon(coord, True, linewidth=0.))
+        fesom_temp_hr.append(selm.var)
+    fesom_salt_hr = []
+    for selm in selements_salt_hr:
+        fesom_salt_hr.append(selm.var)
+    lat_reg = linspace(lat_min, lat_max, num_lat)
+    temp_reg_hr = zeros([num_depth, num_lat])
+    salt_reg_hr = zeros([num_depth, num_lat])
+    temp_reg_hr[:,:] = NaN
+    salt_reg_hr[:,:] = NaN
+    for elm in elements_hr:
+        if amin(elm.lon) < lon0 and amax(elm.lon) > lon0:
+            if amax(elm.lat) > lat_min and amin(elm.lat) < lat_max:
+                tmp = nonzero(lat_reg > amin(elm.lat))[0]
+                if len(tmp) == 0:
+                    jS = 0
+                else:
+                    jS = tmp[0] - 1
+                tmp = nonzero(lat_reg > amax(elm.lat))[0]
+                if len(tmp) == 0:
+                    jN = num_lat
+                else:
+                    jN = tmp[0]
+                for j in range(jS+1,jN):
+                    lat0 = lat_reg[j]
+                    if in_triangle(elm, lon0, lat0):
+                        area = triangle_area(elm.lon, elm.lat)
+                        area0 = triangle_area([lon0, elm.lon[1], elm.lon[2]], [lat0, elm.lat[1], elm.lat[2]])
+                        area1 = triangle_area([lon0, elm.lon[0], elm.lon[2]], [lat0, elm.lat[0], elm.lat[2]])
+                        area2 = triangle_area([lon0, elm.lon[0], elm.lon[1]], [lat0, elm.lat[0], elm.lat[1]])
+                        cff = [area0/area, area1/area, area2/area]
+                        for k in range(num_depth):
+                            node_vals_temp = []
+                            node_vals_salt = []
+                            for n in range(3):
+                                id1, id2, coeff1, coeff2 = elm.nodes[n].find_depth(depth_reg[k])
+                                if any(isnan(array([id1, id2, coeff1, coeff2]))):
+                                    node_vals_temp.append(NaN)
+                                    node_vals_salt.append(NaN)
+                                else:
+                                    node_vals_temp.append(coeff1*fesom_temp_nodes_hr[id1] + coeff2*fesom_temp_nodes_hr[id2])
+                                    node_vals_salt.append(coeff1*fesom_salt_nodes_hr[id1] + coeff2*fesom_salt_nodes_hr[id2])
+                            if any(isnan(node_vals_temp)):
+                                pass
+                            else:
+                                temp_reg_hr[k,j] = sum(array(cff)*array(node_vals_temp))
+                                salt_reg_hr[k,j] = sum(array(cff)*array(node_vals_salt))
+    temp_reg_hr = ma.masked_where(isnan(temp_reg_hr), temp_reg_hr)
+    salt_reg_hr = ma.masked_where(isnan(salt_reg_hr), salt_reg_hr)
+
     depth_reg = -1*depth_reg
 
     # Set up axis labels the way we want them
@@ -209,10 +274,10 @@ def mip_drift_slices (roms_grid, roms_file, fesom_mesh_path, fesom_file):
         depth_labels.append(str(int(round(-val))))
 
     print 'Plotting'
-    fig = figure(figsize=(14,18))
+    fig = figure(figsize=(14,24))
     # ECCO2
     gs1 = GridSpec(1,2)
-    gs1.update(left=0.1, right=0.95, bottom=0.69, top=0.93, wspace=0.08)
+    gs1.update(left=0.1, right=0.95, bottom=0.7575, top=0.94, wspace=0.08)
     # Temperature
     ax = subplot(gs1[0,0])
     pcolor(ecco_lat, ecco_depth, ecco_temp, vmin=temp_min, vmax=temp_max, cmap='jet')
@@ -240,7 +305,7 @@ def mip_drift_slices (roms_grid, roms_file, fesom_mesh_path, fesom_file):
     ax.set_yticklabels([])
     # MetROMS
     gs2 = GridSpec(1,2)
-    gs2.update(left=0.1, right=0.95, bottom=0.38, top=0.62, wspace=0.08)
+    gs2.update(left=0.1, right=0.95, bottom=0.525, top=0.7075, wspace=0.08)
     # Temperature
     ax = subplot(gs2[0,0])
     pcolor(roms_lat, roms_z, roms_temp, vmin=temp_min, vmax=temp_max, cmap='jet')
@@ -263,18 +328,51 @@ def mip_drift_slices (roms_grid, roms_file, fesom_mesh_path, fesom_file):
     ax.set_xticklabels(lat_labels, fontsize=16)
     ax.set_yticks(depth_ticks)
     ax.set_yticklabels([])
-    # FESOM
+    # FESOM low-res
     gs3 = GridSpec(1,2)
-    gs3.update(left=0.1, right=0.95, bottom=0.07, top=0.31, wspace=0.08)
+    gs3.update(left=0.1, right=0.95, bottom=0.2925, top=0.475, wspace=0.08)
     # Temperature
     ax = subplot(gs3[0,0])
-    img = PatchCollection(patches, cmap='jet')
-    img.set_array(array(fesom_temp))
+    img = PatchCollection(patches_lr, cmap='jet')
+    img.set_array(array(fesom_temp_lr))
     img.set_edgecolor('face')
     img.set_clim(vmin=temp_min, vmax=temp_max)
     ax.add_collection(img)
     # Overlay contour on regular grid
-    contour(lat_reg, depth_reg, temp_reg, levels=[temp_contour], color='black')
+    contour(lat_reg, depth_reg, temp_reg_lr, levels=[temp_contour], color='black')
+    ylabel('Depth (m)', fontsize=18)
+    xlim([lat_min, lat_max])
+    ylim([depth_min, depth_max])
+    ax.set_xticks(lat_ticks)
+    ax.set_xticklabels(lat_labels, fontsize=16)
+    ax.set_yticks(depth_ticks)
+    ax.set_yticklabels(depth_labels, fontsize=16)
+    text(-53, 300, 'c) FESOM (low-res), January 2016', fontsize=28)
+    # Salinity
+    ax = subplot(gs3[0,1])
+    img = PatchCollection(patches_lr, cmap='jet')
+    img.set_array(array(fesom_salt_lr))
+    img.set_edgecolor('face')
+    img.set_clim(vmin=salt_min, vmax=salt_max)
+    ax.add_collection(img)
+    contour(lat_reg, depth_reg, salt_reg_lr, levels=[salt_contour], color='black')
+    xlim([lat_min, lat_max])
+    ylim([depth_min, depth_max])
+    ax.set_xticks(lat_ticks)
+    ax.set_xticklabels(lat_labels, fontsize=16)
+    ax.set_yticks(depth_ticks)
+    ax.set_yticklabels([])
+    # FESOM high-res
+    gs4 = GridSpec(1,2)
+    gs4.update(left=0.1, right=0.95, bottom=0.06, top=0.2425, wspace=0.08)
+    # Temperature
+    ax = subplot(gs4[0,0])
+    img = PatchCollection(patches_hr, cmap='jet')
+    img.set_array(array(fesom_temp_hr))
+    img.set_edgecolor('face')
+    img.set_clim(vmin=temp_min, vmax=temp_max)
+    ax.add_collection(img)
+    contour(lat_reg, depth_reg, temp_reg_hr, levels=[temp_contour], color='black')
     ylabel('Depth (m)', fontsize=18)
     xlim([lat_min, lat_max])
     ylim([depth_min, depth_max])
@@ -282,19 +380,19 @@ def mip_drift_slices (roms_grid, roms_file, fesom_mesh_path, fesom_file):
     ax.set_xticklabels(lat_labels, fontsize=16)
     ax.set_yticks(depth_ticks)
     ax.set_yticklabels(depth_labels, fontsize=16)
-    text(-53, 300, 'c) FESOM (high-res), January 2016', fontsize=28)
+    text(-53, 300, 'd) FESOM (high-res), January 2016', fontsize=28)
     # Add a colorbar for temperature
-    cbaxes = fig.add_axes([0.17, 0.02, 0.3, 0.02])
+    cbaxes = fig.add_axes([0.17, 0.015, 0.3, 0.015])
     cbar = colorbar(img, orientation='horizontal', cax=cbaxes, extend='both', ticks=arange(temp_min, temp_max+2, 2))
     cbar.ax.tick_params(labelsize=16)
     # Salinity
-    ax = subplot(gs3[0,1])
-    img = PatchCollection(patches, cmap='jet')
-    img.set_array(array(fesom_salt))
+    ax = subplot(gs4[0,1])
+    img = PatchCollection(patches_hr, cmap='jet')
+    img.set_array(array(fesom_salt_hr))
     img.set_edgecolor('face')
     img.set_clim(vmin=salt_min, vmax=salt_max)
     ax.add_collection(img)
-    contour(lat_reg, depth_reg, salt_reg, levels=[salt_contour], color='black')
+    contour(lat_reg, depth_reg, salt_reg_hr, levels=[salt_contour], color='black')
     xlim([lat_min, lat_max])
     ylim([depth_min, depth_max])
     ax.set_xticks(lat_ticks)
@@ -302,7 +400,7 @@ def mip_drift_slices (roms_grid, roms_file, fesom_mesh_path, fesom_file):
     ax.set_yticks(depth_ticks)
     ax.set_yticklabels([])
     # Add a colorbar for salinity
-    cbaxes = fig.add_axes([0.6, 0.02, 0.3, 0.02])
+    cbaxes = fig.add_axes([0.6, 0.015, 0.3, 0.02])
     cbar = colorbar(img, orientation='horizontal', cax=cbaxes, extend='both', ticks=arange(salt_min+0.1, salt_max+0.1, 0.2))
     cbar.ax.tick_params(labelsize=16)
     fig.show()
@@ -314,8 +412,10 @@ def mip_drift_slices (roms_grid, roms_file, fesom_mesh_path, fesom_file):
 
     roms_grid = raw_input("Path to ROMS grid file: ")
     roms_file = raw_input("Path to ROMS file containing monthly averaged temperature and salinity for January 2016: ")
-    fesom_mesh_path = raw_input("Path to FESOM mesh directory: ")
-    fesom_file = raw_input("Path to FESOM file containing monthly averaged temperature and salinity for January 2016: ")
-    mip_drift_slices(roms_grid, roms_file, fesom_mesh_path, fesom_file)
+    fesom_mesh_path_lr = raw_input("Path to FESOM low-res mesh directory: ")
+    fesom_file_lr = raw_input("Path to FESOM low-res file containing monthly averaged temperature and salinity for January 2016: ")
+    fesom_mesh_path_hr = raw_input("Path to FESOM high-res mesh directory: ")
+    fesom_file_hr = raw_input("Path to FESOM high-res file containing monthly averaged temperature and salinity for January 2016: ")
+    mip_drift_slices(roms_grid, roms_file, fesom_mesh_path_lr, fesom_file_lr, fesom_mesh_path_hr, fesom_file_hr)
     
     
diff --git a/mip_hi_seasonal.py b/mip_hi_seasonal.py
index 350deb7..50164b4 100644
--- a/mip_hi_seasonal.py
+++ b/mip_hi_seasonal.py
@@ -12,11 +12,13 @@
 # Input:
 # cice_seasonal_file = path to seasonal climatology of CICE variables 'aice' and
 #                      'hi', pre-computed using seasonal_climatology_cice.py
-# fesom_mesh_path = path to FESOM mesh directories
-# fesom_seasonal_file = path to seasonal climatology of FESOM variable 'hice', 
-#                       pre-computed using seasonal_climatology.py in the
-#                       "fesomtools" repository
-def mip_hi_seasonal (cice_seasonal_file, fesom_mesh_path, fesom_seasonal_file):
+# fesom_mesh_path_lr, fesom_mesh_path_hr = paths to FESOM mesh directories for
+#                     low-res and high-res mesh respectively
+# fesom_seasonal_file_lr, fesom_seasonal_file_hr = paths to seasonal
+#                         climatologies of FESOM variable 'hice', pre-computed
+#                         using seasonal_climatology.py in the "fesomtools"
+#                         repository, for low-res and high-res respectively
+def mip_hi_seasonal (cice_seasonal_file, fesom_mesh_path_lr, fesom_seasonal_file_lr, fesom_mesh_path_hr, fesom_seasonal_file_hr):
 
     # Boundaries on plot (under polar coordinate transformation)
     x_min = -36.25
@@ -61,33 +63,49 @@ def mip_hi_seasonal (cice_seasonal_file, fesom_mesh_path, fesom_seasonal_file):
     cice_x = -(cice_lat+90)*cos(cice_lon*deg2rad+pi/2)
     cice_y = (cice_lat+90)*sin(cice_lon*deg2rad+pi/2)
 
-    print 'Processing FESOM'
+    print 'Processing low-res FESOM'
     # Build FESOM mesh
-    elements, patches = make_patches(fesom_mesh_path, circumpolar, mask_cavities)
+    elements_lr, patches_lr = make_patches(fesom_mesh_path_lr, circumpolar, mask_cavities)
     # Read data
-    id = Dataset(fesom_seasonal_file, 'r')
-    fesom_hi_nodes = id.variables['hice'][:,:]
+    id = Dataset(fesom_seasonal_file_lr, 'r')
+    fesom_hi_nodes_lr = id.variables['hice'][:,:]
     id.close()
     # Count the number of elements not in ice shelf cavities
-    num_elm = 0
-    for elm in elements:
+    num_elm_lr = 0
+    for elm in elements_lr:
         if not elm.cavity:
-            num_elm += 1
+            num_elm_lr += 1
     # Set up array for element-averages for each season
-    fesom_hi = zeros([4, num_elm])
+    fesom_hi_lr = zeros([4, num_elm_lr])
     # Loop over elements to fill this in
     i = 0
-    for elm in elements:
+    for elm in elements_lr:
         if not elm.cavity:
             # Average over 3 component nodes
-            fesom_hi[:,i] = (fesom_hi_nodes[:,elm.nodes[0].id] + fesom_hi_nodes[:,elm.nodes[1].id] + fesom_hi_nodes[:,elm.nodes[2].id])/3
+            fesom_hi_lr[:,i] = (fesom_hi_nodes_lr[:,elm.nodes[0].id] + fesom_hi_nodes_lr[:,elm.nodes[1].id] + fesom_hi_nodes_lr[:,elm.nodes[2].id])/3
+            i += 1
+
+    print 'Processing high-res FESOM'
+    elements_hr, patches_hr = make_patches(fesom_mesh_path_hr, circumpolar, mask_cavities)
+    id = Dataset(fesom_seasonal_file_hr, 'r')
+    fesom_hi_nodes_hr = id.variables['hice'][:,:]
+    id.close()
+    num_elm_hr = 0
+    for elm in elements_hr:
+        if not elm.cavity:
+            num_elm_hr += 1
+    fesom_hi_hr = zeros([4, num_elm_hr])
+    i = 0
+    for elm in elements_hr:
+        if not elm.cavity:
+            fesom_hi_hr[:,i] = (fesom_hi_nodes_hr[:,elm.nodes[0].id] + fesom_hi_nodes_hr[:,elm.nodes[1].id] + fesom_hi_nodes_hr[:,elm.nodes[2].id])/3
             i += 1
 
     print 'Plotting'
-    fig = figure(figsize=(19,9))
+    fig = figure(figsize=(19,14))
     for season in range(4):
         # MetROMS
-        ax = fig.add_subplot(2, 4, season+1, aspect='equal')
+        ax = fig.add_subplot(3, 4, season+1, aspect='equal')
         pcolor(cice_x, cice_y, cice_hi[season,:,:], vmin=bounds[0], vmax=bounds[1], cmap='jet')
         if season == 0:
             text(-43, 0, 'MetROMS', fontsize=24, ha='right')
@@ -96,10 +114,24 @@ def mip_hi_seasonal (cice_seasonal_file, fesom_mesh_path, fesom_seasonal_file):
         ylim([y_min, y_max])
         ax.set_xticks([])
         ax.set_yticks([])
-        # FESOM
-        ax = fig.add_subplot(2, 4, season+5, aspect='equal')
-        img = PatchCollection(patches, cmap='jet')
-        img.set_array(fesom_hi[season,:])
+        # FESOM low-res
+        ax = fig.add_subplot(3, 4, season+5, aspect='equal')
+        img = PatchCollection(patches_lr, cmap='jet')
+        img.set_array(fesom_hi_lr[season,:])
+        img.set_clim(vmin=bounds[0], vmax=bounds[1])
+        img.set_edgecolor('face')
+        ax.add_collection(img)
+        xlim([x_min, x_max])
+        ylim([y_min, y_max])
+        ax.set_xticks([])
+        ax.set_yticks([])
+        if season == 0:
+            text(-43, 0, 'FESOM', fontsize=24, ha='right')
+            text(-43, -10, '(low-res)', fontsize=24,ha='right')
+        # FESOM high-res
+        ax = fig.add_subplot(3, 4, season+9, aspect='equal')
+        img = PatchCollection(patches_hr, cmap='jet')
+        img.set_array(fesom_hi_hr[season,:])
         img.set_clim(vmin=bounds[0], vmax=bounds[1])
         img.set_edgecolor('face')
         ax.add_collection(img)
@@ -110,11 +142,12 @@ def mip_hi_seasonal (cice_seasonal_file, fesom_mesh_path, fesom_seasonal_file):
         if season == 0:
             text(-43, 0, 'FESOM', fontsize=24, ha='right')
             text(-43, -10, '(high-res)', fontsize=24,ha='right')
-    cbaxes = fig.add_axes([0.35, 0.04, 0.3, 0.02])
+    cbaxes = fig.add_axes([0.35, 0.03, 0.3, 0.02])
     cbar = colorbar(img, orientation='horizontal', ticks=arange(bounds[0],bounds[1]+0.5,0.5), cax=cbaxes, extend='max')
     cbar.ax.tick_params(labelsize=20)
     suptitle('Sea ice effective thickness (m), 1992-2016 average', fontsize=30)
     subplots_adjust(wspace=0.025,hspace=0.025)
+    fig.show()
     fig.savefig('hi_seasonal.png')
 
 
diff --git a/mip_iceshelf_figures.py b/mip_iceshelf_figures.py
index af22d58..105e6eb 100644
--- a/mip_iceshelf_figures.py
+++ b/mip_iceshelf_figures.py
@@ -2,7 +2,7 @@
 from numpy import *
 from numpy.ma import MaskedArray
 from matplotlib.pyplot import *
-from matplotlib.patches import Polygon
+from matplotlib.patches import Polygon, Rectangle
 from matplotlib.collections import PatchCollection, LineCollection
 from matplotlib.cm import *
 from matplotlib.colors import LinearSegmentedColormap, ListedColormap
@@ -264,7 +264,7 @@ def make_vectors (x_min, x_max, y_min, y_max, num_bins_x, num_bins_y):
     # Loop over all points (can't find a better way to do this)
     for j in range(size(roms_speed,0)):
         for i in range(size(roms_speed,1)):
-            # Make sure data isn't masked (i.e. land or open ocean)
+            # Make sure data isn't masked (i.e. land)
             if u_circ_roms[j,i] is not ma.masked:
                 # Check if we're in the region of interest
                 if roms_x[j,i] > x_min and roms_x[j,i] < x_max and roms_y[j,i] > y_min and roms_y[j,i] < y_max:
@@ -293,13 +293,12 @@ def make_vectors (x_min, x_max, y_min, y_max, num_bins_x, num_bins_y):
     v_circ_fesom_lr = node_speed_lr*sin(theta_circ_fesom_lr)
     # Loop over 2D nodes to fill in the velocity bins as before
     for n in range(fesom_n2d_lr):
-        if fesom_cavity_lr[n]:
-            if fesom_x_lr[n] > x_min and fesom_x_lr[n] < x_max and fesom_y_lr[n] > y_min and fesom_y_lr[n] < y_max:
-                x_index = nonzero(x_bins > fesom_x_lr[n])[0][0]-1
-                y_index = nonzero(y_bins > fesom_y_lr[n])[0][0]-1
-                fesom_ubin_lr[y_index, x_index] += u_circ_fesom_lr[n]
-                fesom_vbin_lr[y_index, x_index] += v_circ_fesom_lr[n]
-                fesom_num_pts_lr[y_index, x_index] += 1
+        if fesom_x_lr[n] > x_min and fesom_x_lr[n] < x_max and fesom_y_lr[n] > y_min and fesom_y_lr[n] < y_max:
+            x_index = nonzero(x_bins > fesom_x_lr[n])[0][0]-1
+            y_index = nonzero(y_bins > fesom_y_lr[n])[0][0]-1
+            fesom_ubin_lr[y_index, x_index] += u_circ_fesom_lr[n]
+            fesom_vbin_lr[y_index, x_index] += v_circ_fesom_lr[n]
+            fesom_num_pts_lr[y_index, x_index] += 1
     fesom_ubin_lr = ma.masked_where(fesom_num_pts_lr==0, fesom_ubin_lr)
     fesom_vbin_lr = ma.masked_where(fesom_num_pts_lr==0, fesom_vbin_lr)
     flag = fesom_num_pts_lr > 0
@@ -314,13 +313,12 @@ def make_vectors (x_min, x_max, y_min, y_max, num_bins_x, num_bins_y):
     u_circ_fesom_hr = node_speed_hr*cos(theta_circ_fesom_hr)
     v_circ_fesom_hr = node_speed_hr*sin(theta_circ_fesom_hr)
     for n in range(fesom_n2d_hr):
-        if fesom_cavity_hr[n]:
-            if fesom_x_hr[n] > x_min and fesom_x_hr[n] < x_max and fesom_y_hr[n] > y_min and fesom_y_hr[n] < y_max:
-                x_index = nonzero(x_bins > fesom_x_hr[n])[0][0]-1
-                y_index = nonzero(y_bins > fesom_y_hr[n])[0][0]-1
-                fesom_ubin_hr[y_index, x_index] += u_circ_fesom_hr[n]
-                fesom_vbin_hr[y_index, x_index] += v_circ_fesom_hr[n]
-                fesom_num_pts_hr[y_index, x_index] += 1
+        if fesom_x_hr[n] > x_min and fesom_x_hr[n] < x_max and fesom_y_hr[n] > y_min and fesom_y_hr[n] < y_max:
+            x_index = nonzero(x_bins > fesom_x_hr[n])[0][0]-1
+            y_index = nonzero(y_bins > fesom_y_hr[n])[0][0]-1
+            fesom_ubin_hr[y_index, x_index] += u_circ_fesom_hr[n]
+            fesom_vbin_hr[y_index, x_index] += v_circ_fesom_hr[n]
+            fesom_num_pts_hr[y_index, x_index] += 1
     fesom_ubin_hr = ma.masked_where(fesom_num_pts_hr==0, fesom_ubin_hr)
     fesom_vbin_hr = ma.masked_where(fesom_num_pts_hr==0, fesom_vbin_hr)
     flag = fesom_num_pts_hr > 0
@@ -495,7 +493,9 @@ def plot_wct (x_min, x_max, y_min, y_max, gs, cbaxes, cbar_ticks, letter):
 # lon_lines = list of longitudes to overlay as dotted lines on the first panel
 #             (-180 to 180)
 # lat_lines = list of latitudes to overlay (-90 to 90)
-def plot_melt (x_min, x_max, y_min, y_max, gs, cbaxes, cbar_ticks, change_points, letter, y0, lon_lines, lat_lines):
+# rect_bounds = optional list of size 4 containing x_min, x_max, y_min, y_max
+#               of a rectangle to superimpose on the third panel
+def plot_melt (x_min, x_max, y_min, y_max, gs, cbaxes, cbar_ticks, change_points, letter, y0, lon_lines, lat_lines, rect_bounds=None):
 
     # Set up a grey square for FESOM to fill the background with land
     x_reg_fesom, y_reg_fesom = meshgrid(linspace(x_min, x_max, num=100), linspace(y_min, y_max, num=100))
@@ -586,10 +586,13 @@ def plot_melt (x_min, x_max, y_min, y_max, gs, cbaxes, cbar_ticks, change_points
     ylim([y_min, y_max])
     ax.set_xticks([])
     ax.set_yticks([])
+    if rect_bounds is not None:
+        # Overlay rectangle
+        ax.add_patch(Rectangle((rect_bounds[0], rect_bounds[2]), rect_bounds[1]-rect_bounds[0], rect_bounds[3]-rect_bounds[2], fill=False))
     # Colourbar
     cbar = colorbar(img, cax=cbaxes, ticks=cbar_ticks)
     text(0.5, y0, 'FESOM (high-res)', fontsize=18, horizontalalignment='center', transform=ax.transAxes) 
-
+    
 
 # Plot bottom water temperature for MetROMS, low-res FESOM, and high-res FESOM
 # for the given region, in the given axes.
@@ -785,7 +788,7 @@ def plot_velavg (x_min, x_max, y_min, y_max, gs, cbaxes, cbar_ticks, x_centres,
     x_reg_fesom, y_reg_fesom = meshgrid(linspace(x_min, x_max, num=100), linspace(y_min, y_max, num=100))
     land_square = zeros(shape(x_reg_fesom))
     # Find bounds on variable in this region
-    var_min, var_max = get_min_max(roms_speed, fesom_speed_lr, fesom_speed_hr, x_min, x_max, y_min, y_max)
+    var_min, var_max = get_min_max(roms_speed, fesom_speed_lr, fesom_speed_hr, x_min, x_max, y_min, y_max, cavity=False)
 
     # Plot MetROMS
     ax = subplot(gs[0,0], aspect='equal')
@@ -796,9 +799,12 @@ def plot_velavg (x_min, x_max, y_min, y_max, gs, cbaxes, cbar_ticks, x_centres,
     pcolor(x_reg_roms, y_reg_roms, land_circle, cmap=grey_cmap)
     # Now shade the data
     pcolor(roms_x, roms_y, roms_speed, vmin=var_min, vmax=var_max, cmap='cool')
+    # Add a black contour for the ice shelf front
+    rcParams['contour.negative_linestyle'] = 'solid'
     if x_centres is not None:
         # Overlay vectors    
         quiver(x_centres, y_centres, roms_ubin, roms_vbin, scale=arrow_scale, headwidth=arrow_headwidth, headlength=arrow_headlength, color='black')
+    contour(roms_x, roms_y, zice_contour, levels=[min_zice], colors=('black'))
     xlim([x_min, x_max])
     ylim([y_min, y_max])
     ax.set_xticks([])
@@ -811,15 +817,14 @@ def plot_velavg (x_min, x_max, y_min, y_max, gs, cbaxes, cbar_ticks, x_centres,
     # Start with land background
     contourf(x_reg_fesom, y_reg_fesom, land_square, 1, colors=(('0.6', '0.6', '0.6')))
     # Add ice shelf elements
-    img = PatchCollection(patches_lr, cmap='cool')
+    img = PatchCollection(patches_all_lr, cmap='cool')
     img.set_array(array(fesom_speed_lr))
     img.set_edgecolor('face')
     img.set_clim(vmin=var_min, vmax=var_max)
     ax.add_collection(img)
-    # Mask out the open ocean in white
-    overlay = PatchCollection(mask_patches_lr, facecolor=(1,1,1))
-    overlay.set_edgecolor('face')
-    ax.add_collection(overlay)
+    # Add ice shelf front contour lines
+    fesom_contours_lr = LineCollection(contour_lines_lr, edgecolor='black', linewidth=1)
+    ax.add_collection(fesom_contours_lr)
     if x_centres is not None:
         # Overlay vectors
         quiver(x_centres, y_centres, fesom_ubin_lr, fesom_vbin_lr, scale=arrow_scale, headwidth=arrow_headwidth, headlength=arrow_headlength, color='black')
@@ -838,14 +843,13 @@ def plot_velavg (x_min, x_max, y_min, y_max, gs, cbaxes, cbar_ticks, x_centres,
     # FESOM high-res
     ax = subplot(gs[0,2], aspect='equal')
     contourf(x_reg_fesom, y_reg_fesom, land_square, 1, colors=(('0.6', '0.6', '0.6')))
-    img = PatchCollection(patches_hr, cmap='cool')
+    img = PatchCollection(patches_all_hr, cmap='cool')
     img.set_array(array(fesom_speed_hr))
     img.set_edgecolor('face')
     img.set_clim(vmin=var_min, vmax=var_max)
     ax.add_collection(img)
-    overlay = PatchCollection(mask_patches_hr, facecolor=(1,1,1))
-    overlay.set_edgecolor('face')
-    ax.add_collection(overlay)
+    fesom_contours_hr = LineCollection(contour_lines_hr, edgecolor='black', linewidth=1)
+    ax.add_collection(fesom_contours_hr)
     if x_centres is not None:
         quiver(x_centres, y_centres, fesom_ubin_hr, fesom_vbin_hr, scale=arrow_scale, headwidth=arrow_headwidth, headlength=arrow_headlength, color='black')
     xlim([x_min, x_max])
@@ -1047,8 +1051,8 @@ def plot_zonal_ts (lon0, lat_min, lat_max, lat_ticks, lat_labels, depth_min, dep
 # File paths
 roms_grid = '/short/m68/kaa561/metroms_iceshelf/apps/common/grid/circ30S_quarterdegree.nc'
 roms_file = '/short/m68/kaa561/metroms_iceshelf/tmproms/run/intercomparison/2002_2016_avg.nc'
-fesom_mesh_path_lr = '/short/y99/kaa561/FESOM/mesh/low_res/'
-fesom_mesh_path_hr = '/short/y99/kaa561/FESOM/mesh/high_res/'
+fesom_mesh_path_lr = '/short/y99/kaa561/FESOM/mesh/meshA/'
+fesom_mesh_path_hr = '/short/y99/kaa561/FESOM/mesh/meshB/'
 fesom_file_lr_o = '/short/y99/kaa561/FESOM/intercomparison_lowres/output/oce_2002_2016_avg.nc'
 fesom_file_hr_o = '/short/y99/kaa561/FESOM/intercomparison_highres/output/oce_2002_2016_avg.nc'
 fesom_file_lr_i = '/short/y99/kaa561/FESOM/intercomparison_lowres/output/wnet_2002_2016_avg.nc'
@@ -1327,9 +1331,9 @@ def plot_zonal_ts (lon0, lat_min, lat_max, lat_ticks, lat_labels, depth_min, dep
 # Vertically average u and v
 roms_u = sum(u_3d*dz, axis=0)/sum(dz, axis=0)
 roms_v = sum(v_3d*dz, axis=0)/sum(dz, axis=0)
-# Mask the open ocean and land
-roms_u = ma.masked_where(roms_zice==0, roms_u)
-roms_v = ma.masked_where(roms_zice==0, roms_v)
+# Apply land mask
+roms_u = ma.masked_where(roms_mask==0, roms_u)
+roms_v = ma.masked_where(roms_mask==0, roms_v)
 # Calculate speed
 roms_speed = sqrt(roms_u**2 + roms_v**2)
 
@@ -1424,8 +1428,7 @@ def plot_zonal_ts (lon0, lat_min, lat_max, lat_ticks, lat_labels, depth_min, dep
 # Calculate speed at each element, averaged over 3 corners
 fesom_speed_lr = []
 for elm in elements_lr:
-    if elm.cavity:
-        fesom_speed_lr.append(mean([node_speed_lr[elm.nodes[0].id], node_speed_lr[elm.nodes[1].id], node_speed_lr[elm.nodes[2].id]]))
+    fesom_speed_lr.append(mean([node_speed_lr[elm.nodes[0].id], node_speed_lr[elm.nodes[1].id], node_speed_lr[elm.nodes[2].id]]))
 
 # FESOM high-res
 fesom_cavity_hr = []
@@ -1497,8 +1500,7 @@ def plot_zonal_ts (lon0, lat_min, lat_max, lat_ticks, lat_labels, depth_min, dep
 node_speed_hr = sqrt(node_u_hr**2 + node_v_hr**2)
 fesom_speed_hr = []
 for elm in elements_hr:
-    if elm.cavity:
-        fesom_speed_hr.append(mean([node_speed_hr[elm.nodes[0].id], node_speed_hr[elm.nodes[1].id], node_speed_hr[elm.nodes[2].id]])) 
+    fesom_speed_hr.append(mean([node_speed_hr[elm.nodes[0].id], node_speed_hr[elm.nodes[1].id], node_speed_hr[elm.nodes[2].id]])) 
 
 
 # **************** USER MODIFIED SECTION ****************
@@ -1507,45 +1509,39 @@ def plot_zonal_ts (lon0, lat_min, lat_max, lat_ticks, lat_labels, depth_min, dep
 x_max_tmp = -4.5
 y_min_tmp = 1
 y_max_tmp = 10
-fig = figure(figsize=(8,22))
+fig = figure(figsize=(8,14))
 fig.patch.set_facecolor('white')
 # Melt
 gs_a = GridSpec(1,3)
-gs_a.update(left=0.05, right=0.9, bottom=0.764, top=0.89, wspace=0.05)
-cbaxes_tmp = fig.add_axes([0.91, 0.777, 0.025, 0.1])
+gs_a.update(left=0.05, right=0.9, bottom=0.735, top=0.89, wspace=0.05)
+cbaxes_tmp = fig.add_axes([0.91, 0.755, 0.025, 0.12])
 cbar_ticks = arange(0, 6+3, 3)
 plot_melt(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_a, cbaxes_tmp, cbar_ticks, [0.5, 3, 4.5], 'a', 1.25, [-80, -60, -40], [-80])
 # Velocity
 x_centres, y_centres, roms_ubin, roms_vbin, fesom_ubin_lr, fesom_vbin_lr, fesom_ubin_hr, fesom_vbin_hr = make_vectors(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, 20, 20)
 gs_b = GridSpec(1,3)
-gs_b.update(left=0.05, right=0.9, bottom=0.614, top=0.74, wspace=0.05)
-cbaxes_tmp = fig.add_axes([0.91, 0.627, 0.025, 0.11])
+gs_b.update(left=0.05, right=0.9, bottom=0.555, top=0.71, wspace=0.05)
+cbaxes_tmp = fig.add_axes([0.91, 0.575, 0.025, 0.12])
 cbar_ticks = arange(0, 0.2+0.1, 0.1)
 plot_velavg(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_b, cbaxes_tmp, cbar_ticks, x_centres, y_centres, roms_ubin, roms_vbin, fesom_ubin_lr, fesom_vbin_lr, fesom_ubin_hr, fesom_vbin_hr, 'b', arrow_scale=0.9, arrow_headwidth=8, arrow_headlength=9)
 # Bottom water temperature
 gs_c = GridSpec(1,3)
-gs_c.update(left=0.05, right=0.9, bottom=0.464, top=0.59, wspace=0.05)
-cbaxes_tmp = fig.add_axes([0.91, 0.477, 0.025, 0.1])
-cbar_ticks = arange(-2.6, -1.0+0.8, 0.8)
+gs_c.update(left=0.05, right=0.9, bottom=0.375, top=0.53, wspace=0.05)
+cbaxes_tmp = fig.add_axes([0.91, 0.395, 0.025, 0.12])
+cbar_ticks = arange(-2.4, -1.0+0.7, 0.7)
 plot_bwtemp(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_c, cbaxes_tmp, cbar_ticks, 'c')
 # Bottom water salinity
 gs_d = GridSpec(1,3)
-gs_d.update(left=0.05, right=0.9, bottom=0.314, top=0.44, wspace=0.05)
-cbaxes_tmp = fig.add_axes([0.91, 0.327, 0.025, 0.1])
+gs_d.update(left=0.05, right=0.9, bottom=0.195, top=0.35, wspace=0.05)
+cbaxes_tmp = fig.add_axes([0.91, 0.215, 0.025, 0.12])
 cbar_ticks = arange(34.3, 34.7+0.2, 0.2)
 plot_bwsalt(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_d, cbaxes_tmp, cbar_ticks, 'd')
-# Water column thickness
-gs_e = GridSpec(1,3)
-gs_e.update(left=0.05, right=0.9, bottom=0.164, top=0.29, wspace=0.05)
-cbaxes_tmp = fig.add_axes([0.91, 0.177, 0.025, 0.1])
-cbar_ticks = arange(200, 1000+400, 400)
-plot_wct(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_e, cbaxes_tmp, cbar_ticks, 'e')
 # Ice shelf draft
-gs_f = GridSpec(1,3)
-gs_f.update(left=0.05, right=0.9, bottom=0.014, top=0.14, wspace=0.05)
-cbaxes_tmp = fig.add_axes([0.91, 0.027, 0.025, 0.1])
+gs_e = GridSpec(1,3)
+gs_e.update(left=0.05, right=0.9, bottom=0.015, top=0.17, wspace=0.05)
+cbaxes_tmp = fig.add_axes([0.91, 0.035, 0.025, 0.12])
 cbar_ticks = arange(500, 1500+500, 500)
-plot_draft(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_f, cbaxes_tmp, cbar_ticks, 'f')
+plot_draft(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_e, cbaxes_tmp, cbar_ticks, 'e')
 suptitle('Filchner-Ronne Ice Shelf', fontsize=30)
 fig.show()
 fig.savefig('filchner_ronne.png')
@@ -1564,12 +1560,12 @@ def plot_zonal_ts (lon0, lat_min, lat_max, lat_ticks, lat_labels, depth_min, dep
 cbar_ticks = arange(0, 4+2, 2)
 plot_melt(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_a, cbaxes_tmp, cbar_ticks, [1, 2, 3], 'a', 1.35, [-10, 20], [-70])
 # Velocity
-#x_centres, y_centres, roms_ubin, roms_vbin, fesom_ubin_lr, fesom_vbin_lr, fesom_ubin_hr, fesom_vbin_hr = make_vectors(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, 40, 20)
+x_centres, y_centres, roms_ubin, roms_vbin, fesom_ubin_lr, fesom_vbin_lr, fesom_ubin_hr, fesom_vbin_hr = make_vectors(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, 20, 10)
 gs_b = GridSpec(1,3)
 gs_b.update(left=0.11, right=0.9, bottom=0.57, top=0.72, wspace=0.05)
 cbaxes_tmp = fig.add_axes([0.91, 0.6, 0.025, 0.1])
 cbar_ticks = arange(0, 0.18+0.09, 0.09)
-plot_velavg(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_b, cbaxes_tmp, cbar_ticks, None, None, None, None, None, None, None, None, 'b', arrow_scale=0.9, arrow_headwidth=8, arrow_headlength=9)
+plot_velavg(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_b, cbaxes_tmp, cbar_ticks, x_centres, y_centres, roms_ubin, roms_vbin, fesom_ubin_lr, fesom_vbin_lr, fesom_ubin_hr, fesom_vbin_hr, 'b', arrow_scale=0.9, arrow_headwidth=8, arrow_headlength=9)
 # Temp and salt slices through 1W: Fimbul Ice Shelf
 lat_min = -71.5
 lat_max = -69.5
@@ -1597,45 +1593,39 @@ def plot_zonal_ts (lon0, lat_min, lat_max, lat_ticks, lat_labels, depth_min, dep
 x_max_tmp = 20.5
 y_min_tmp = 4.75
 y_max_tmp = 8
-fig = figure(figsize=(10,14))
+fig = figure(figsize=(10,12))
 fig.patch.set_facecolor('white')
 # Melt
 gs_a = GridSpec(1,3)
-gs_a.update(left=0.05, right=0.9, bottom=0.764, top=0.89, wspace=0.05)
-cbaxes_tmp = fig.add_axes([0.91, 0.777, 0.025, 0.1])
+gs_a.update(left=0.05, right=0.9, bottom=0.735, top=0.89, wspace=0.05)
+cbaxes_tmp = fig.add_axes([0.91, 0.755, 0.025, 0.12])
 cbar_ticks = arange(0, 40+20, 20)
 plot_melt(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_a, cbaxes_tmp, cbar_ticks, [2, 5, 11], 'a', 1.25, [70], [-70])
 # Velocity
 x_centres, y_centres, roms_ubin, roms_vbin, fesom_ubin_lr, fesom_vbin_lr, fesom_ubin_hr, fesom_vbin_hr = make_vectors(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, 20, 15)
 gs_b = GridSpec(1,3)
-gs_b.update(left=0.05, right=0.9, bottom=0.614, top=0.74, wspace=0.05)
-cbaxes_tmp = fig.add_axes([0.91, 0.627, 0.025, 0.11])
+gs_b.update(left=0.05, right=0.9, bottom=0.555, top=0.71, wspace=0.05)
+cbaxes_tmp = fig.add_axes([0.91, 0.575, 0.025, 0.12])
 cbar_ticks = arange(0.1, 0.3+0.1, 0.1)
 plot_velavg(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_b, cbaxes_tmp, cbar_ticks, x_centres, y_centres, roms_ubin, roms_vbin, fesom_ubin_lr, fesom_vbin_lr, fesom_ubin_hr, fesom_vbin_hr, 'b', arrow_scale=0.9, arrow_headwidth=7, arrow_headlength=8)
 # Bottom water temperature
 gs_c = GridSpec(1,3)
-gs_c.update(left=0.05, right=0.9, bottom=0.464, top=0.59, wspace=0.05)
-cbaxes_tmp = fig.add_axes([0.91, 0.477, 0.025, 0.1])
+gs_c.update(left=0.05, right=0.9, bottom=0.375, top=0.53, wspace=0.05)
+cbaxes_tmp = fig.add_axes([0.91, 0.395, 0.025, 0.12])
 cbar_ticks = arange(-2.4, -1.4+0.5, 0.5)
 plot_bwtemp(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_c, cbaxes_tmp, cbar_ticks, 'c')
 # Bottom water salinity
 gs_d = GridSpec(1,3)
-gs_d.update(left=0.05, right=0.9, bottom=0.314, top=0.44, wspace=0.05)
-cbaxes_tmp = fig.add_axes([0.91, 0.327, 0.025, 0.1])
-cbar_ticks = arange(34, 34.4+0.2, 0.2)
+gs_d.update(left=0.05, right=0.9, bottom=0.195, top=0.35, wspace=0.05)
+cbaxes_tmp = fig.add_axes([0.91, 0.215, 0.025, 0.12])
+cbar_ticks = arange(34.2, 34.4+0.2, 0.2)
 plot_bwsalt(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_d, cbaxes_tmp, cbar_ticks, 'd')
-# Water column thickness
-gs_e = GridSpec(1,3)
-gs_e.update(left=0.05, right=0.9, bottom=0.164, top=0.29, wspace=0.05)
-cbaxes_tmp = fig.add_axes([0.91, 0.177, 0.025, 0.1])
-cbar_ticks = arange(250, 1250+500, 500)
-plot_wct(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_e, cbaxes_tmp, cbar_ticks, 'e')
 # Ice shelf draft
-gs_f = GridSpec(1,3)
-gs_f.update(left=0.05, right=0.9, bottom=0.014, top=0.14, wspace=0.05)
-cbaxes_tmp = fig.add_axes([0.91, 0.027, 0.025, 0.1])
+gs_e = GridSpec(1,3)
+gs_e.update(left=0.05, right=0.9, bottom=0.015, top=0.17, wspace=0.05)
+cbaxes_tmp = fig.add_axes([0.91, 0.035, 0.025, 0.12])
 cbar_ticks = arange(200, 2200+1000, 1000)
-plot_draft(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_f, cbaxes_tmp, cbar_ticks, 'f')
+plot_draft(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_e, cbaxes_tmp, cbar_ticks, 'e')
 suptitle('Amery Ice Shelf', fontsize=30)
 fig.show()
 fig.savefig('amery.png')
@@ -1645,6 +1635,11 @@ def plot_zonal_ts (lon0, lat_min, lat_max, lat_ticks, lat_labels, depth_min, dep
 x_max_tmp = 25.5
 y_min_tmp = -20
 y_max_tmp = 4
+# Zoomed-in region for Totten
+x_min_totten = 20
+x_max_totten = 21.5
+y_min_totten = -10.9
+y_max_totten = -9.7
 fig = figure(figsize=(8,14))
 fig.patch.set_facecolor('white')
 # Melt
@@ -1652,31 +1647,28 @@ def plot_zonal_ts (lon0, lat_min, lat_max, lat_ticks, lat_labels, depth_min, dep
 gs_a.update(left=0.05, right=0.9, bottom=0.68, top=0.9, wspace=0.05)
 cbaxes_tmp = fig.add_axes([0.91, 0.72, 0.025, 0.15])
 cbar_ticks = arange(0, 4+2, 2)
-plot_melt(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_a, cbaxes_tmp, cbar_ticks, [1, 2, 3], 'a', 1.15, [90, 120], [-65])
+# Overlay rectangle showing Totten region to zoom in
+plot_melt(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_a, cbaxes_tmp, cbar_ticks, [1, 2, 3], 'a', 1.15, [90, 120], [-65], rect_bounds=[x_min_totten, x_max_totten, y_min_totten, y_max_totten])
 # Bottom water temperature
 x_reg, y_reg, fesom_bathy_reg_lr, fesom_bathy_reg_hr = fesom_interpolate_bathy(75, 160, -70, -63)
 gs_b = GridSpec(1,3)
 gs_b.update(left=0.05, right=0.9, bottom=0.435, top=0.655, wspace=0.05)
 cbaxes_tmp = fig.add_axes([0.91, 0.475, 0.025, 0.15])
 cbar_ticks = arange(-1.5, 0.5+1, 1)
-plot_bwtemp(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_b, cbaxes_tmp, cbar_ticks, 'b', bathy_contour=1000)
+plot_bwtemp(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_b, cbaxes_tmp, cbar_ticks, 'b', bathy_contour=1500)
 # Bottom water salinity
 gs_c = GridSpec(1,3)
 gs_c.update(left=0.05, right=0.9, bottom=0.19, top=0.41, wspace=0.05)
 cbaxes_tmp = fig.add_axes([0.91, 0.23, 0.025, 0.15])
 cbar_ticks = arange(34.2, 34.6+0.2, 0.2)
-plot_bwsalt(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_c, cbaxes_tmp, cbar_ticks, 'c', bathy_contour=1000)
+plot_bwsalt(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_c, cbaxes_tmp, cbar_ticks, 'c', bathy_contour=1500)
 # Velocity for Totten
-x_min_tmp = 20
-x_max_tmp = 21.5
-y_min_tmp = -10.9
-y_max_tmp = -9.7
-x_centres, y_centres, roms_ubin, roms_vbin, fesom_ubin_lr, fesom_vbin_lr, fesom_ubin_hr, fesom_vbin_hr = make_vectors(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, 18, 18)
+x_centres, y_centres, roms_ubin, roms_vbin, fesom_ubin_lr, fesom_vbin_lr, fesom_ubin_hr, fesom_vbin_hr = make_vectors(x_min_totten, x_max_totten, y_min_totten, y_max_totten, 18, 18)
 gs_d = GridSpec(1,3)
 gs_d.update(left=0.05, right=0.9, bottom=0.025, top=0.165, wspace=0.05)
 cbaxes_tmp = fig.add_axes([0.91, 0.05, 0.025, 0.08])
-cbar_ticks = arange(0.02, 0.08+0.03, 0.03)
-plot_velavg(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_d, cbaxes_tmp, cbar_ticks, x_centres, y_centres, roms_ubin, roms_vbin, fesom_ubin_lr, fesom_vbin_lr, fesom_ubin_hr, fesom_vbin_hr, 'd', loc_string='Totten Ice Shelf', arrow_scale=0.7, arrow_headwidth=8, arrow_headlength=9)
+cbar_ticks = arange(0.01, 0.07+0.03, 0.03)
+plot_velavg(x_min_totten, x_max_totten, y_min_totten, y_max_totten, gs_d, cbaxes_tmp, cbar_ticks, x_centres, y_centres, roms_ubin, roms_vbin, fesom_ubin_lr, fesom_vbin_lr, fesom_ubin_hr, fesom_vbin_hr, 'd', loc_string='Totten Ice Shelf', arrow_scale=0.7, arrow_headwidth=8, arrow_headlength=9)
 suptitle('Australian Sector', fontsize=30)
 fig.show()
 fig.savefig('australian.png')
@@ -1686,27 +1678,21 @@ def plot_zonal_ts (lon0, lat_min, lat_max, lat_ticks, lat_labels, depth_min, dep
 x_max_tmp = 4
 y_min_tmp = -13
 y_max_tmp = -4.75
-fig = figure(figsize=(10,12))
+fig = figure(figsize=(10,10))
 fig.patch.set_facecolor('white')
 # Melt
 gs_a = GridSpec(1,3)
-gs_a.update(left=0.1, right=0.9, bottom=0.735, top=0.89, wspace=0.05)
-cbaxes_tmp = fig.add_axes([0.91, 0.765, 0.025, 0.1])
+gs_a.update(left=0.1, right=0.9, bottom=0.69, top=0.9, wspace=0.05)
+cbaxes_tmp = fig.add_axes([0.91, 0.73, 0.025, 0.13])
 cbar_ticks = arange(0, 8+4, 4)
 plot_melt(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_a, cbaxes_tmp, cbar_ticks, [0.5, 2, 4], 'a', 1.25, [180, -140], [-80])
 # Velocity
 x_centres, y_centres, roms_ubin, roms_vbin, fesom_ubin_lr, fesom_vbin_lr, fesom_ubin_hr, fesom_vbin_hr = make_vectors(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, 20, 15)
 gs_b = GridSpec(1,3)
-gs_b.update(left=0.1, right=0.9, bottom=0.56, top=0.715, wspace=0.05)
-cbaxes_tmp = fig.add_axes([0.91, 0.59, 0.025, 0.1])
+gs_b.update(left=0.1, right=0.9, bottom=0.47, top=0.68, wspace=0.05)
+cbaxes_tmp = fig.add_axes([0.91, 0.51, 0.025, 0.13])
 cbar_ticks = arange(0, 0.18+0.09, 0.09)
 plot_velavg(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_b, cbaxes_tmp, cbar_ticks, x_centres, y_centres, roms_ubin, roms_vbin, fesom_ubin_lr, fesom_vbin_lr, fesom_ubin_hr, fesom_vbin_hr, 'b', arrow_scale=0.5, arrow_headwidth=7, arrow_headlength=8)
-# Draft
-gs_c = GridSpec(1,3)
-gs_c.update(left=0.1, right=0.9, bottom=0.385, top=0.54, wspace=0.05)
-cbaxes_tmp = fig.add_axes([0.91, 0.415, 0.025, 0.1])
-cbar_ticks = arange(200, 600+200, 200)
-plot_draft(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_c, cbaxes_tmp, cbar_ticks, 'c')
 # Temp and salt slices through 180E
 lat_min = -84.5
 lat_max = -77
@@ -1716,24 +1702,29 @@ def plot_zonal_ts (lon0, lat_min, lat_max, lat_ticks, lat_labels, depth_min, dep
 depth_max = 0
 depth_ticks = [-750, -500, -250, 0]
 depth_labels = ['750', '500', '250', '0']
-gs_d = GridSpec(1,3)
-gs_d.update(left=0.1, right=0.9, bottom=0.23, top=0.36, wspace=0.05)
-cbaxes_tmp1 = fig.add_axes([0.91, 0.245, 0.025, 0.1])
+gs_c = GridSpec(1,3)
+gs_c.update(left=0.1, right=0.9, bottom=0.28, top=0.44, wspace=0.05)
+cbaxes_tmp1 = fig.add_axes([0.91, 0.3, 0.025, 0.12])
 cbar_ticks1 = arange(-2.1, -1.5+0.3, 0.3)
-gs_e = GridSpec(1,3)
-gs_e.update(left=0.1, right=0.9, bottom=0.05, top=0.18, wspace=0.05)
-cbaxes_tmp2 = fig.add_axes([0.91, 0.065, 0.025, 0.1])
+gs_d = GridSpec(1,3)
+gs_d.update(left=0.1, right=0.9, bottom=0.05, top=0.21, wspace=0.05)
+cbaxes_tmp2 = fig.add_axes([0.91, 0.07, 0.025, 0.12])
 cbar_ticks2 = arange(34.4, 35+0.3, 0.3)
-plot_zonal_ts(180, lat_min, lat_max, lat_ticks, lat_labels, depth_min, depth_max, depth_ticks, depth_labels, gs_d, gs_e, cbaxes_tmp1, cbaxes_tmp2, cbar_ticks1, cbar_ticks2, 'd', 'e')
+plot_zonal_ts(180, lat_min, lat_max, lat_ticks, lat_labels, depth_min, depth_max, depth_ticks, depth_labels, gs_c, gs_d, cbaxes_tmp1, cbaxes_tmp2, cbar_ticks1, cbar_ticks2, 'c', 'd')
 suptitle('Ross Sea', fontsize=30)
 fig.show()
 fig.savefig('ross.png')
 
 # Amundsen Sea
-x_min_tmp = -15.5
+x_min_tmp = -17.5
 x_max_tmp = -10.5
 y_min_tmp = -11.25
 y_max_tmp = -2.25
+# Zoomed-in region for PIG
+x_min_pig = -15.6
+x_max_pig = -14.1
+y_min_pig = -3.5
+y_max_pig = -2.25
 fig = figure(figsize=(8,14))
 fig.patch.set_facecolor('white')
 # Melt
@@ -1741,31 +1732,27 @@ def plot_zonal_ts (lon0, lat_min, lat_max, lat_ticks, lat_labels, depth_min, dep
 gs_a.update(left=0.05, right=0.9, bottom=0.68, top=0.9, wspace=0.05)
 cbaxes_tmp = fig.add_axes([0.91, 0.72, 0.025, 0.15])
 cbar_ticks = arange(0, 14+7, 7)
-plot_melt(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_a, cbaxes_tmp, cbar_ticks, [1, 3, 6], 'a', 1.15, [-130, -110], [-75])
+plot_melt(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_a, cbaxes_tmp, cbar_ticks, [1, 3, 6], 'a', 1.15, [-130, -110], [-75], rect_bounds=[x_min_pig, x_max_pig, y_min_pig, y_max_pig])
 # Bottom water temperature
-x_reg, y_reg, fesom_bathy_reg_lr, fesom_bathy_reg_hr = fesom_interpolate_bathy(-140, -90, -76, -72)
+x_reg, y_reg, fesom_bathy_reg_lr, fesom_bathy_reg_hr = fesom_interpolate_bathy(-140, -90, -76, -70)
 gs_b = GridSpec(1,3)
 gs_b.update(left=0.05, right=0.9, bottom=0.435, top=0.655, wspace=0.05)
 cbaxes_tmp = fig.add_axes([0.91, 0.475, 0.025, 0.15])
 cbar_ticks = arange(-1.6, 0.8+0.8, 0.8)
-plot_bwtemp(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_b, cbaxes_tmp, cbar_ticks, 'b', bathy_contour=1000)
+plot_bwtemp(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_b, cbaxes_tmp, cbar_ticks, 'b', bathy_contour=1500)
 # Bottom water salinity
 gs_c = GridSpec(1,3)
 gs_c.update(left=0.05, right=0.9, bottom=0.19, top=0.41, wspace=0.05)
 cbaxes_tmp = fig.add_axes([0.91, 0.23, 0.025, 0.15])
 cbar_ticks = arange(34.1, 34.7+0.2, 0.2)
-plot_bwsalt(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_c, cbaxes_tmp, cbar_ticks, 'c', bathy_contour=1000)
+plot_bwsalt(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_c, cbaxes_tmp, cbar_ticks, 'c', bathy_contour=1500)
 # Velocity for PIG
-x_min_tmp = -15.6
-x_max_tmp = -14.1
-y_min_tmp = -3.5
-y_max_tmp = -2.25
-x_centres, y_centres, roms_ubin, roms_vbin, fesom_ubin_lr, fesom_vbin_lr, fesom_ubin_hr, fesom_vbin_hr = make_vectors(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, 18, 18)
+x_centres, y_centres, roms_ubin, roms_vbin, fesom_ubin_lr, fesom_vbin_lr, fesom_ubin_hr, fesom_vbin_hr = make_vectors(x_min_pig, x_max_pig, y_min_pig, y_max_pig, 18, 18)
 gs_d = GridSpec(1,3)
 gs_d.update(left=0.05, right=0.9, bottom=0.025, top=0.165, wspace=0.05)
 cbaxes_tmp = fig.add_axes([0.91, 0.05, 0.025, 0.08])
 cbar_ticks = arange(0.03, 0.09+0.03, 0.03)
-plot_velavg(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_d, cbaxes_tmp, cbar_ticks, x_centres, y_centres, roms_ubin, roms_vbin, fesom_ubin_lr, fesom_vbin_lr, fesom_ubin_hr, fesom_vbin_hr, 'd', loc_string='Pine Island Glacier Ice Shelf', arrow_scale=0.5, arrow_headwidth=8, arrow_headlength=9)
+plot_velavg(x_min_pig, x_max_pig, y_min_pig, y_max_pig, gs_d, cbaxes_tmp, cbar_ticks, x_centres, y_centres, roms_ubin, roms_vbin, fesom_ubin_lr, fesom_vbin_lr, fesom_ubin_hr, fesom_vbin_hr, 'd', loc_string='Pine Island Glacier Ice Shelf', arrow_scale=0.5, arrow_headwidth=8, arrow_headlength=9)
 suptitle('Amundsen Sea', fontsize=30)
 fig.show()
 fig.savefig('amundsen.png')
@@ -1775,6 +1762,11 @@ def plot_zonal_ts (lon0, lat_min, lat_max, lat_ticks, lat_labels, depth_min, dep
 x_max_tmp = -15.5
 y_min_tmp = -4.5
 y_max_tmp = 7.6
+# Zoomed-in region for George VI
+x_min_gvi = -18.75
+x_max_gvi = -15.5
+y_min_gvi = 4.25
+y_max_gvi = 7.6
 fig = figure(figsize=(8,14))
 fig.patch.set_facecolor('white')
 # Melt
@@ -1782,31 +1774,27 @@ def plot_zonal_ts (lon0, lat_min, lat_max, lat_ticks, lat_labels, depth_min, dep
 gs_a.update(left=0.05, right=0.9, bottom=0.68, top=0.9, wspace=0.05)
 cbaxes_tmp = fig.add_axes([0.91, 0.72, 0.025, 0.15])
 cbar_ticks = arange(0, 12+6, 6)
-plot_melt(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_a, cbaxes_tmp, cbar_ticks, [0.5, 2, 4], 'a', 1.15, [-100, -80], [-70])
+plot_melt(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_a, cbaxes_tmp, cbar_ticks, [0.5, 2, 4], 'a', 1.15, [-100, -80], [-70], rect_bounds=[x_min_gvi, x_max_gvi, y_min_gvi, y_max_gvi])
 # Bottom water temperature
 x_reg, y_reg, fesom_bathy_reg_lr, fesom_bathy_reg_hr = fesom_interpolate_bathy(-110, -60, -76, -66)
 gs_b = GridSpec(1,3)
 gs_b.update(left=0.05, right=0.9, bottom=0.435, top=0.655, wspace=0.05)
 cbaxes_tmp = fig.add_axes([0.91, 0.475, 0.025, 0.15])
 cbar_ticks = arange(-1.5, 0.5+1, 1)
-plot_bwtemp(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_b, cbaxes_tmp, cbar_ticks, 'b', bathy_contour=1000)
+plot_bwtemp(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_b, cbaxes_tmp, cbar_ticks, 'b', bathy_contour=1500)
 # Bottom water salinity
 gs_c = GridSpec(1,3)
 gs_c.update(left=0.05, right=0.9, bottom=0.19, top=0.41, wspace=0.05)
 cbaxes_tmp = fig.add_axes([0.91, 0.23, 0.025, 0.15])
 cbar_ticks = arange(33.8, 34.6+0.4, 0.4)
-plot_bwsalt(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_c, cbaxes_tmp, cbar_ticks, 'c', bathy_contour=1000)
+plot_bwsalt(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_c, cbaxes_tmp, cbar_ticks, 'c', bathy_contour=1500)
 # Velocity for George VI
-x_min_tmp = -18.75
-x_max_tmp = -15.5
-y_min_tmp = 4.25
-y_max_tmp = 7.6
-x_centres, y_centres, roms_ubin, roms_vbin, fesom_ubin_lr, fesom_vbin_lr, fesom_ubin_hr, fesom_vbin_hr = make_vectors(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, 18, 18)
+x_centres, y_centres, roms_ubin, roms_vbin, fesom_ubin_lr, fesom_vbin_lr, fesom_ubin_hr, fesom_vbin_hr = make_vectors(x_min_gvi, x_max_gvi, y_min_gvi, y_max_gvi, 18, 18)
 gs_d = GridSpec(1,3)
 gs_d.update(left=0.05, right=0.9, bottom=0.025, top=0.165, wspace=0.05)
 cbaxes_tmp = fig.add_axes([0.91, 0.05, 0.025, 0.08])
-cbar_ticks = arange(0, 0.08+0.04, 0.04)
-plot_velavg(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_d, cbaxes_tmp, cbar_ticks, x_centres, y_centres, roms_ubin, roms_vbin, fesom_ubin_lr, fesom_vbin_lr, fesom_ubin_hr, fesom_vbin_hr, 'd', loc_string='George VI Ice Shelf', arrow_scale=0.4, arrow_headwidth=8, arrow_headlength=9)
+cbar_ticks = arange(0, 0.1+0.05, 0.05)
+plot_velavg(x_min_gvi, x_max_gvi, y_min_gvi, y_max_gvi, gs_d, cbaxes_tmp, cbar_ticks, x_centres, y_centres, roms_ubin, roms_vbin, fesom_ubin_lr, fesom_vbin_lr, fesom_ubin_hr, fesom_vbin_hr, 'd', loc_string='George VI Ice Shelf', arrow_scale=0.4, arrow_headwidth=8, arrow_headlength=9)
 suptitle('Bellingshausen Sea', fontsize=30)
 fig.show()
 fig.savefig('bellingshausen.png')
@@ -1822,21 +1810,22 @@ def plot_zonal_ts (lon0, lat_min, lat_max, lat_ticks, lat_labels, depth_min, dep
 gs_a = GridSpec(1,3)
 gs_a.update(left=0.05, right=0.9, bottom=0.61, top=0.84, wspace=0.05)
 cbaxes_tmp = fig.add_axes([0.91, 0.65, 0.025, 0.15])
-cbar_ticks = arange(0, 8+4, 4)
+cbar_ticks = arange(0, 6+3, 3)
 plot_melt(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_a, cbaxes_tmp, cbar_ticks, [0.5, 2, 4], 'a', 1.3, [-60], [-70])
 # Velocity
 x_centres, y_centres, roms_ubin, roms_vbin, fesom_ubin_lr, fesom_vbin_lr, fesom_ubin_hr, fesom_vbin_hr = make_vectors(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, 18, 18)
 gs_b = GridSpec(1,3)
 gs_b.update(left=0.05, right=0.9, bottom=0.33, top=0.56, wspace=0.05)
 cbaxes_tmp = fig.add_axes([0.91, 0.37, 0.025, 0.15])
-cbar_ticks = arange(0, 0.12+0.06, 0.06)
+cbar_ticks = arange(0, 0.2+0.1, 0.1)
 plot_velavg(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_b, cbaxes_tmp, cbar_ticks, x_centres, y_centres, roms_ubin, roms_vbin, fesom_ubin_lr, fesom_vbin_lr, fesom_ubin_hr, fesom_vbin_hr, 'b', arrow_scale=0.7, arrow_headwidth=8, arrow_headlength=9)
 # Bottom water temperature
+x_reg, y_reg, fesom_bathy_reg_lr, fesom_bathy_reg_hr = fesom_interpolate_bathy(-70, -40, -75, -63)
 gs_c = GridSpec(1,3)
 gs_c.update(left=0.05, right=0.9, bottom=0.05, top=0.28, wspace=0.05)
 cbaxes_tmp = fig.add_axes([0.91, 0.09, 0.025, 0.15])
 cbar_ticks = arange(-1.8, 0+0.9, 0.9)
-plot_bwtemp(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_c, cbaxes_tmp, cbar_ticks, 'c')
+plot_bwtemp(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_c, cbaxes_tmp, cbar_ticks, 'c', bathy_contour=1500)
 suptitle('Larsen Ice Shelves', fontsize=30)
 fig.show()
 fig.savefig('larsen.png')
diff --git a/mip_mld.py b/mip_mld.py
index ced2e21..0c03da9 100644
--- a/mip_mld.py
+++ b/mip_mld.py
@@ -15,11 +15,14 @@
 # roms_grid = path to ROMS grid file
 # roms_seasonal_file = path to seasonal climatology of ROMS 3D temperature and
 #                      salinity, precomputed using seasonal_climatology_roms.py
-# fesom_mesh_path = path to FESOM mesh directory
-# fesom_seasonal_file = path to seasonal climatology of FESOM 3D temperature and
-#                       salinity, precomputed using seasonal_climatology.py in
-#                       the "fesomtools" repository
-def mip_mld (roms_grid, roms_seasonal_file, fesom_mesh_path, fesom_seasonal_file):
+# fesom_mesh_path_lr, fesom_mesh_path_hr = path to FESOM mesh directories for
+#                     low-res and high-res meshes
+# fesom_seasonal_file_lr, fesom_seasonal_file_hr = paths to seasonal
+#                         climatologies of FESOM 3D temperature and salinity
+#                         for low-res and high-res respectively, precomputed
+#                         using seasonal_climatology.py in the "fesomtools"
+#                         repository
+def mip_mld (roms_grid, roms_seasonal_file, fesom_mesh_path_lr, fesom_seasonal_file_lr, fesom_mesh_path_hr, fesom_seasonal_file_hr):
 
     # Path to Sallee's observations
     obs_file = '/short/m68/kaa561/Climatology_MLD003_v2017.nc'
@@ -106,29 +109,29 @@ def mip_mld (roms_grid, roms_seasonal_file, fesom_mesh_path, fesom_seasonal_file
                             break
                         k -= 1
 
-    print 'Processing FESOM:'
+    print 'Processing low-res FESOM:'
     print 'Building mesh'
-    elements, patches = make_patches(fesom_mesh_path, circumpolar, mask_cavities)
+    elements_lr, patches_lr = make_patches(fesom_mesh_path_lr, circumpolar, mask_cavities)
     print 'Reading data'
-    id = Dataset(fesom_seasonal_file, 'r')
-    fesom_temp_nodes = id.variables['temp'][:,:]
-    fesom_salt_nodes = id.variables['salt'][:,:]
+    id = Dataset(fesom_seasonal_file_lr, 'r')
+    fesom_temp_nodes_lr = id.variables['temp'][:,:]
+    fesom_salt_nodes_lr = id.variables['salt'][:,:]
     id.close()
     print 'Calculating density'
-    fesom_density_nodes = unesco(fesom_temp_nodes, fesom_salt_nodes, zeros(shape(fesom_temp_nodes)))
+    fesom_density_nodes_lr = unesco(fesom_temp_nodes_lr, fesom_salt_nodes_lr, zeros(shape(fesom_temp_nodes_lr)))
     print 'Calculating mixed layer depth'
     # Set up array for mixed layer depth at each element, at each season
-    fesom_mld = zeros([4, len(elements)])
+    fesom_mld_lr = zeros([4, len(elements_lr)])
     # Loop over seasons and elements to fill these in
     for season in range(4):
         print '...' + season_names[season]
         mld_season = []
-        for elm in elements:
+        for elm in elements_lr:
             # Get mixed layer depth at each node
             mld_nodes = []
             for i in range(3):
                 node = elm.nodes[i]
-                density_sfc = fesom_density_nodes[season,node.id]
+                density_sfc = fesom_density_nodes_lr[season,node.id]
                 # Save surface depth (only nonzero in ice shelf cavities)
                 depth_sfc = node.depth
                 temp_depth = node.depth
@@ -138,7 +141,7 @@ def mip_mld (roms_grid, roms_seasonal_file, fesom_mesh_path, fesom_seasonal_file
                         # Reached the bottom
                         mld_nodes.append(temp_depth-depth_sfc)
                         break
-                    if fesom_density_nodes[season,curr_node.id] >= density_sfc + density_anom:
+                    if fesom_density_nodes_lr[season,curr_node.id] >= density_sfc + density_anom:
                         # Reached the critical density anomaly
                         mld_nodes.append(curr_node.depth-depth_sfc)
                         break
@@ -146,7 +149,49 @@ def mip_mld (roms_grid, roms_seasonal_file, fesom_mesh_path, fesom_seasonal_file
                     curr_node = curr_node.below
             # For this element, save the mean mixed layer depth
             mld_season.append(mean(array(mld_nodes)))
-        fesom_mld[season,:] = array(mld_season)
+        fesom_mld_lr[season,:] = array(mld_season)
+
+    print 'Processing high-res FESOM:'
+    print 'Building mesh'
+    elements_hr, patches_hr = make_patches(fesom_mesh_path_hr, circumpolar, mask_cavities)
+    print 'Reading data'
+    id = Dataset(fesom_seasonal_file_hr, 'r')
+    fesom_temp_nodes_hr = id.variables['temp'][:,:]
+    fesom_salt_nodes_hr = id.variables['salt'][:,:]
+    id.close()
+    print 'Calculating density'
+    fesom_density_nodes_hr = unesco(fesom_temp_nodes_hr, fesom_salt_nodes_hr, zeros(shape(fesom_temp_nodes_hr)))
+    print 'Calculating mixed layer depth'
+    # Set up array for mixed layer depth at each element, at each season
+    fesom_mld_hr = zeros([4, len(elements_hr)])
+    # Loop over seasons and elements to fill these in
+    for season in range(4):
+        print '...' + season_names[season]
+        mld_season = []
+        for elm in elements_hr:
+            # Get mixed layer depth at each node
+            mld_nodes = []
+            for i in range(3):
+                node = elm.nodes[i]
+                density_sfc = fesom_density_nodes_hr[season,node.id]
+                # Save surface depth (only nonzero in ice shelf cavities)
+                depth_sfc = node.depth
+                temp_depth = node.depth
+                curr_node = node.below
+                while True:
+                    if curr_node is None:
+                        # Reached the bottom
+                        mld_nodes.append(temp_depth-depth_sfc)
+                        break
+                    if fesom_density_nodes_hr[season,curr_node.id] >= density_sfc + density_anom:
+                        # Reached the critical density anomaly
+                        mld_nodes.append(curr_node.depth-depth_sfc)
+                        break
+                    temp_depth = curr_node.depth
+                    curr_node = curr_node.below
+            # For this element, save the mean mixed layer depth
+            mld_season.append(mean(array(mld_nodes)))
+        fesom_mld_hr[season,:] = array(mld_season)
 
     print 'Processing obs'
     # Read grid and monthly climatology
@@ -185,12 +230,12 @@ def mip_mld (roms_grid, roms_seasonal_file, fesom_mesh_path, fesom_seasonal_file
 
     print 'Plotting'
     # ACC
-    fig1 = figure(figsize=(13,9))
+    fig1 = figure(figsize=(18,9))
     # Summer
     # MetROMS
-    ax = fig1.add_subplot(2, 3, 1, aspect='equal')
+    ax = fig1.add_subplot(2, 4, 1, aspect='equal')
     pcolor(roms_x, roms_y, roms_mld[0,:,:], vmin=0, vmax=max_bound_summer, cmap='jet')
-    text(-65, 0, season_names[0], fontsize=24, ha='right')
+    text(-67, 0, season_names[0], fontsize=24, ha='right')
     title('MetROMS', fontsize=24)
     xlim([-nbdry1, nbdry1])
     ylim([-nbdry1, nbdry1])
@@ -199,10 +244,22 @@ def mip_mld (roms_grid, roms_seasonal_file, fesom_mesh_path, fesom_seasonal_file
         text(x_ticks[i], y_ticks[i], lon_labels[i], ha='center', rotation=lon_rot[i], fontsize=12)
     ax.set_xticks([])
     ax.set_yticks([])
-    # FESOM
-    ax = fig1.add_subplot(2, 3, 2, aspect='equal')
-    img = PatchCollection(patches, cmap='jet')
-    img.set_array(fesom_mld[0,:])
+    # FESOM low-res
+    ax = fig1.add_subplot(2, 4, 2, aspect='equal')
+    img = PatchCollection(patches_lr, cmap='jet')
+    img.set_array(fesom_mld_lr[0,:])
+    img.set_clim(vmin=0, vmax=max_bound_summer)
+    img.set_edgecolor('face')
+    ax.add_collection(img)
+    xlim([-nbdry1, nbdry1])
+    ylim([-nbdry1, nbdry1])
+    ax.set_xticks([])
+    ax.set_yticks([])
+    title('FESOM (low-res)', fontsize=24)
+    # FESOM high-res
+    ax = fig1.add_subplot(2, 4, 3, aspect='equal')
+    img = PatchCollection(patches_hr, cmap='jet')
+    img.set_array(fesom_mld_hr[0,:])
     img.set_clim(vmin=0, vmax=max_bound_summer)
     img.set_edgecolor('face')
     ax.add_collection(img)
@@ -212,7 +269,7 @@ def mip_mld (roms_grid, roms_seasonal_file, fesom_mesh_path, fesom_seasonal_file
     ax.set_yticks([])
     title('FESOM (high-res)', fontsize=24)
     # Obs
-    ax = fig1.add_subplot(2, 3, 3, aspect='equal')
+    ax = fig1.add_subplot(2, 4, 4, aspect='equal')
     img = pcolor(obs_x, obs_y, obs_mld[0,:,:], vmin=0, vmax=max_bound_summer, cmap='jet')    
     xlim([-nbdry1, nbdry1])
     ylim([-nbdry1, nbdry1])
@@ -225,17 +282,28 @@ def mip_mld (roms_grid, roms_seasonal_file, fesom_mesh_path, fesom_seasonal_file
     cbar.ax.tick_params(labelsize=20)
     # Winter
     # MetROMS
-    ax = fig1.add_subplot(2, 3, 4, aspect='equal')
+    ax = fig1.add_subplot(2, 4, 5, aspect='equal')
     pcolor(roms_x, roms_y, roms_mld[2,:,:], vmin=0, vmax=max_bound_winter, cmap='jet')
-    text(-65, 0, season_names[2], fontsize=24, ha='right')
+    text(-67, 0, season_names[2], fontsize=24, ha='right')
+    xlim([-nbdry1, nbdry1])
+    ylim([-nbdry1, nbdry1])
+    ax.set_xticks([])
+    ax.set_yticks([])
+    # FESOM low-res
+    ax = fig1.add_subplot(2, 4, 6, aspect='equal')
+    img = PatchCollection(patches_lr, cmap='jet')
+    img.set_array(fesom_mld_lr[2,:])
+    img.set_clim(vmin=0, vmax=max_bound_winter)
+    img.set_edgecolor('face')
+    ax.add_collection(img)
     xlim([-nbdry1, nbdry1])
     ylim([-nbdry1, nbdry1])
     ax.set_xticks([])
     ax.set_yticks([])
-    # FESOM
-    ax = fig1.add_subplot(2, 3, 5, aspect='equal')
-    img = PatchCollection(patches, cmap='jet')
-    img.set_array(fesom_mld[2,:])
+    # FESOM high-res
+    ax = fig1.add_subplot(2, 4, 7, aspect='equal')
+    img = PatchCollection(patches_hr, cmap='jet')
+    img.set_array(fesom_mld_hr[2,:])
     img.set_clim(vmin=0, vmax=max_bound_winter)
     img.set_edgecolor('face')
     ax.add_collection(img)
@@ -244,7 +312,7 @@ def mip_mld (roms_grid, roms_seasonal_file, fesom_mesh_path, fesom_seasonal_file
     ax.set_xticks([])
     ax.set_yticks([])
     # Obs
-    ax = fig1.add_subplot(2, 3, 6, aspect='equal')
+    ax = fig1.add_subplot(2, 4, 8, aspect='equal')
     img = pcolor(obs_x, obs_y, obs_mld[2,:,:], vmin=0, vmax=max_bound_winter, cmap='jet')
     xlim([-nbdry1, nbdry1])
     ylim([-nbdry1, nbdry1])
@@ -262,10 +330,10 @@ def mip_mld (roms_grid, roms_seasonal_file, fesom_mesh_path, fesom_seasonal_file
     fig1.savefig('mld_acc.png')
 
     # Continental shelf
-    fig2 = figure(figsize=(9.5,9))
+    fig2 = figure(figsize=(13,9))
     # Summer
     # MetROMS
-    ax = fig2.add_subplot(2, 2, 1, aspect='equal')
+    ax = fig2.add_subplot(2, 3, 1, aspect='equal')
     pcolor(roms_x, roms_y, roms_mld[0,:,:], vmin=0, vmax=max_bound_summer, cmap='jet')
     text(-28, 0, season_names[0], fontsize=24, ha='right')
     title('MetROMS', fontsize=24)
@@ -273,10 +341,22 @@ def mip_mld (roms_grid, roms_seasonal_file, fesom_mesh_path, fesom_seasonal_file
     ylim([-nbdry2, nbdry2])
     ax.set_xticks([])
     ax.set_yticks([])
-    # FESOM
-    ax = fig2.add_subplot(2, 2, 2, aspect='equal')
-    img = PatchCollection(patches, cmap='jet')
-    img.set_array(fesom_mld[0,:])
+    # FESOM low-res
+    ax = fig2.add_subplot(2, 3, 2, aspect='equal')
+    img = PatchCollection(patches_lr, cmap='jet')
+    img.set_array(fesom_mld_lr[0,:])
+    img.set_clim(vmin=0, vmax=max_bound_summer)
+    img.set_edgecolor('face')
+    ax.add_collection(img)
+    xlim([-nbdry2, nbdry2])
+    ylim([-nbdry2, nbdry2])
+    ax.set_xticks([])
+    ax.set_yticks([])
+    title('FESOM (low-res)', fontsize=24)
+    # FESOM high-res
+    ax = fig2.add_subplot(2, 3, 3, aspect='equal')
+    img = PatchCollection(patches_hr, cmap='jet')
+    img.set_array(fesom_mld_hr[0,:])
     img.set_clim(vmin=0, vmax=max_bound_summer)
     img.set_edgecolor('face')
     ax.add_collection(img)
@@ -286,22 +366,33 @@ def mip_mld (roms_grid, roms_seasonal_file, fesom_mesh_path, fesom_seasonal_file
     ax.set_yticks([])
     title('FESOM (high-res)', fontsize=24)
     # Add a colorbar for summer
-    cbaxes = fig2.add_axes([0.91, 0.55, 0.02, 0.3])
+    cbaxes = fig2.add_axes([0.93, 0.55, 0.02, 0.3])
     cbar = colorbar(img, cax=cbaxes, extend='max', ticks=arange(0, max_bound_summer+50, 50))
     cbar.ax.tick_params(labelsize=20)
     # Winter
     # MetROMS
-    ax = fig2.add_subplot(2, 2, 3, aspect='equal')
+    ax = fig2.add_subplot(2, 3, 4, aspect='equal')
     pcolor(roms_x, roms_y, roms_mld[2,:,:], vmin=0, vmax=max_bound_winter, cmap='jet')
     text(-28, 0, season_names[2], fontsize=24, ha='right')
     xlim([-nbdry2, nbdry2])
     ylim([-nbdry2, nbdry2])
     ax.set_xticks([])
     ax.set_yticks([])
-    # FESOM
-    ax = fig2.add_subplot(2, 2, 4, aspect='equal')
-    img = PatchCollection(patches, cmap='jet')
-    img.set_array(fesom_mld[2,:])
+    # FESOM low-res
+    ax = fig2.add_subplot(2, 3, 5, aspect='equal')
+    img = PatchCollection(patches_lr, cmap='jet')
+    img.set_array(fesom_mld_lr[2,:])
+    img.set_clim(vmin=0, vmax=max_bound_winter)
+    img.set_edgecolor('face')
+    ax.add_collection(img)
+    xlim([-nbdry2, nbdry2])
+    ylim([-nbdry2, nbdry2])
+    ax.set_xticks([])
+    ax.set_yticks([])
+    # FESOM high-res
+    ax = fig2.add_subplot(2, 3, 6, aspect='equal')
+    img = PatchCollection(patches_hr, cmap='jet')
+    img.set_array(fesom_mld_hr[2,:])
     img.set_clim(vmin=0, vmax=max_bound_winter)
     img.set_edgecolor('face')
     ax.add_collection(img)
@@ -310,7 +401,7 @@ def mip_mld (roms_grid, roms_seasonal_file, fesom_mesh_path, fesom_seasonal_file
     ax.set_xticks([])
     ax.set_yticks([])
     # Add a colorbar for winter
-    cbaxes = fig2.add_axes([0.91, 0.15, 0.02, 0.3])
+    cbaxes = fig2.add_axes([0.93, 0.15, 0.02, 0.3])
     cbar = colorbar(img, cax=cbaxes, extend='max', ticks=arange(0, max_bound_winter+200, 200))
     cbar.ax.tick_params(labelsize=20)
     # Add the main title
@@ -326,9 +417,11 @@ def mip_mld (roms_grid, roms_seasonal_file, fesom_mesh_path, fesom_seasonal_file
 
     roms_grid = raw_input("Path to ROMS grid file: ")
     roms_seasonal_file = raw_input("Path to ROMS seasonal climatology file containing 3D temp and salt: ")
-    fesom_mesh_path = raw_input("Path to FESOM mesh directory: ")
-    fesom_seasonal_file = raw_input("Path to FESOM seasonal climatology file containing 3D temp and salt: ")
-    mip_mld(roms_grid, roms_seasonal_file, fesom_mesh_path, fesom_seasonal_file)
+    fesom_mesh_path_lr = raw_input("Path to FESOM low-res mesh directory: ")
+    fesom_seasonal_file_lr = raw_input("Path to FESOM low-res seasonal climatology file containing 3D temp and salt: ")
+    fesom_mesh_path_hr = raw_input("Path to FESOM high-res mesh directory: ")
+    fesom_seasonal_file_hr = raw_input("Path to FESOM high-res seasonal climatology file containing 3D temp and salt: ")
+    mip_mld(roms_grid, roms_seasonal_file, fesom_mesh_path_lr, fesom_seasonal_file_lr, fesom_mesh_path_hr, fesom_seasonal_file_hr)
 
     
                             
diff --git a/mip_regions_1var.py b/mip_regions_1var.py
index 4bef5f0..114f3d4 100644
--- a/mip_regions_1var.py
+++ b/mip_regions_1var.py
@@ -11,7 +11,7 @@
 import sys
 sys.path.insert(0, '/short/y99/kaa561/fesomtools')
 from fesom_grid import *
-from patches import *
+from matplotlib.patches import Polygon
 from unrotate_vector import *
 from unrotate_grid import *
 
@@ -30,8 +30,8 @@ def mip_regions_1var ():
     # Path to ROMS time-averaged file
     roms_file = '/short/m68/kaa561/metroms_iceshelf/tmproms/run/intercomparison/2002_2016_avg.nc'
     # Path to FESOM mesh directories
-    fesom_mesh_path_lr = '/short/y99/kaa561/FESOM/mesh/low_res/'
-    fesom_mesh_path_hr = '/short/y99/kaa561/FESOM/mesh/high_res/'
+    fesom_mesh_path_lr = '/short/y99/kaa561/FESOM/mesh/meshA/'
+    fesom_mesh_path_hr = '/short/y99/kaa561/FESOM/mesh/meshB/'
     # Path to FESOM time-averaged ocean files (temp, salt, u, v)
     fesom_file_lr_o = '/short/y99/kaa561/FESOM/intercomparison_lowres/output/oce_2002_2016_avg.nc'
     fesom_file_hr_o = '/short/y99/kaa561/FESOM/intercomparison_highres/output/oce_2002_2016_avg.nc'
@@ -52,7 +52,7 @@ def mip_regions_1var ():
     # Size of each plot in the y direction
     ysize = [8, 6, 7, 9, 7, 9, 10, 7]
     # Variables to process
-    var_names = ['draft', 'melt', 'temp', 'salt', 'vel', 'bathy', 'wct']
+    var_names = ['vel'] #['bathy', 'draft', 'wct', 'melt', 'temp', 'salt', 'vel']
     # Constants
     sec_per_year = 365*24*3600
     deg2rad = pi/180.0
@@ -95,14 +95,31 @@ def mip_regions_1var ():
     land_circle = ma.masked_where(sqrt((x_reg_roms-x_c)**2 + (y_reg_roms-y_c)**2) > radius, land_circle)
 
     print 'Building FESOM low-res mesh'
-    # Mask open ocean
-    elements_lr, mask_patches_lr = make_patches(fesom_mesh_path_lr, circumpolar=True, mask_cavities=True)
-    # Unmask ice shelves
-    patches_lr = iceshelf_mask(elements_lr)
+    elements_lr = fesom_grid(fesom_mesh_path_lr, circumpolar=True)
+    # Make patches for all elements, ice shelf elements, and open ocean elements
+    patches_lr = []
+    patches_shelf_lr = []
+    patches_ocn_lr = []
+    for elm in elements_lr:
+        coord = transpose(vstack((elm.x, elm.y)))
+        patches_lr.append(Polygon(coord, True, linewidth=0.))
+        if elm.cavity:
+            patches_shelf_lr.append(Polygon(coord, True, linewidth=0.))
+        else:
+            patches_ocn_lr.append(Polygon(coord, True, linewidth=0.))
 
     print 'Building FESOM high-res mesh'
-    elements_hr, mask_patches_hr = make_patches(fesom_mesh_path_hr, circumpolar=True, mask_cavities=True)
-    patches_hr = iceshelf_mask(elements_hr)
+    elements_hr = fesom_grid(fesom_mesh_path_hr, circumpolar=True)
+    patches_hr = []
+    patches_shelf_hr = []
+    patches_ocn_hr = []
+    for elm in elements_hr:
+        coord = transpose(vstack((elm.x, elm.y)))
+        patches_hr.append(Polygon(coord, True, linewidth=0.))
+        if elm.cavity:
+            patches_shelf_hr.append(Polygon(coord, True, linewidth=0.))
+        else:
+            patches_ocn_hr.append(Polygon(coord, True, linewidth=0.))
 
     for var in var_names:
         print 'Processing variable ' + var
@@ -112,8 +129,8 @@ def mip_regions_1var ():
             # Swap sign on existing zice field; nothing more to read
             roms_data = -1*roms_zice
         elif var == 'bathy':
-            # Point to h field; nothing more to read
-            roms_data = roms_h
+            # Point to h field and mask out land mask; nothing more to read
+            roms_data = ma.masked_where(roms_mask==0, roms_h)
         elif var == 'wct':
             # Add h (positive) and zice (negative); nothing more to read
             roms_data = roms_h + roms_zice
@@ -170,7 +187,14 @@ def mip_regions_1var ():
                 v_rho = sum(v_3d*dz, axis=0)/sum(dz, axis=0)    
                 # Get speed
                 roms_data = sqrt(u_rho**2 + v_rho**2)
+                # Mask out land
+                u_rho = ma.masked_where(roms_mask==0, u_rho)
+                v_rho = ma.masked_where(roms_mask==0, v_rho)
+                roms_data = ma.masked_where(roms_mask==0, roms_data)
             id.close()
+        if var in ['draft', 'melt', 'wct']:
+            # Mask out open ocean
+            roms_data = ma.masked_where(roms_zice==0, roms_data)
 
         print 'Reading FESOM low-res fields'
         if var not in ['draft', 'bathy', 'wct']:
@@ -274,9 +298,9 @@ def mip_regions_1var ():
         # Calculate given field at each element
         fesom_data_lr = []
         for elm in elements_lr:
-            # For each element in an ice shelf cavity, append the mean value
-            # for the 3 component Nodes
-            if elm.cavity:
+            # For each element, append the mean value for the 3 component Nodes
+            # Restrict to ice shelf cavities for draft, melt, wct
+            if elm.cavity or var not in ['draft', 'melt', 'wct']:
                 if var == 'draft':
                     # Ice shelf draft is depth of surface layer
                     fesom_data_lr.append(mean([elm.nodes[0].depth, elm.nodes[1].depth, elm.nodes[2].depth]))
@@ -374,7 +398,7 @@ def mip_regions_1var ():
             id.close()
         fesom_data_hr = []
         for elm in elements_hr:
-            if elm.cavity:
+            if elm.cavity or var not in ['draft', 'melt', 'wct']:
                 if var == 'draft':
                     fesom_data_hr.append(mean([elm.nodes[0].depth, elm.nodes[1].depth, elm.nodes[2].depth]))
                 elif var == 'bathy':
@@ -401,7 +425,7 @@ def mip_regions_1var ():
             # Low-res
             i = 0
             for elm in elements_lr:
-                if elm.cavity:
+                if elm.cavity or var not in ['draft', 'melt', 'wct']:
                     if any(elm.x >= x_min[index]) and any(elm.x <= x_max[index]) and any(elm.y >= y_min[index]) and any(elm.y <= y_max[index]):
                         if fesom_data_lr[i] < var_min:
                             var_min = fesom_data_lr[i]
@@ -411,7 +435,7 @@ def mip_regions_1var ():
             # High-res
             i = 0
             for elm in elements_hr:
-                if elm.cavity:
+                if elm.cavity or var not in ['draft', 'melt', 'wct']:
                     if any(elm.x >= x_min[index]) and any(elm.x <= x_max[index]) and any(elm.y >= y_min[index]) and any(elm.y <= y_max[index]):
                         if fesom_data_hr[i] < var_min:
                             var_min = fesom_data_hr[i]
@@ -467,7 +491,7 @@ def mip_regions_1var ():
                 # Loop over all points (can't find a better way to do this)
                 for j in range(size(roms_data,0)):
                     for i in range(size(roms_data,1)):
-                        # Make sure data isn't masked (i.e. land or open ocean)
+                        # Make sure data isn't masked (i.e. land)
                         if u_circ_roms[j,i] is not ma.masked:
                             # Check if we're in the region of interest
                             if roms_x[j,i] > x_min[index] and roms_x[j,i] < x_max[index] and roms_y[j,i] > y_min[index] and roms_y[j,i] < y_max[index]:
@@ -496,13 +520,12 @@ def mip_regions_1var ():
                 v_circ_fesom_lr = node_data_lr*sin(theta_circ_fesom_lr)
                 # Loop over 2D nodes to fill in the velocity bins as before
                 for n in range(fesom_n2d_lr):
-                    if fesom_cavity_lr[n]:
-                        if fesom_x_lr[n] > x_min[index] and fesom_x_lr[n] < x_max[index] and fesom_y_lr[n] > y_min[index] and fesom_y_lr[n] < y_max[index]:
-                            x_index = nonzero(x_bins > fesom_x_lr[n])[0][0]-1
-                            y_index = nonzero(y_bins > fesom_y_lr[n])[0][0]-1
-                            fesom_u_lr[y_index, x_index] += u_circ_fesom_lr[n]
-                            fesom_v_lr[y_index, x_index] += v_circ_fesom_lr[n]
-                            fesom_num_pts_lr[y_index, x_index] += 1
+                    if fesom_x_lr[n] > x_min[index] and fesom_x_lr[n] < x_max[index] and fesom_y_lr[n] > y_min[index] and fesom_y_lr[n] < y_max[index]:
+                        x_index = nonzero(x_bins > fesom_x_lr[n])[0][0]-1
+                        y_index = nonzero(y_bins > fesom_y_lr[n])[0][0]-1
+                        fesom_u_lr[y_index, x_index] += u_circ_fesom_lr[n]
+                        fesom_v_lr[y_index, x_index] += v_circ_fesom_lr[n]
+                        fesom_num_pts_lr[y_index, x_index] += 1
                 fesom_u_lr = ma.masked_where(fesom_num_pts_lr==0, fesom_u_lr)
                 fesom_v_lr = ma.masked_where(fesom_num_pts_lr==0, fesom_v_lr)
                 flag = fesom_num_pts_lr > 0
@@ -518,13 +541,12 @@ def mip_regions_1var ():
                 v_circ_fesom_hr = node_data_hr*sin(theta_circ_fesom_hr)
                 # Loop over 2D nodes to fill in the velocity bins as before
                 for n in range(fesom_n2d_hr):
-                    if fesom_cavity_hr[n]:
-                        if fesom_x_hr[n] > x_min[index] and fesom_x_hr[n] < x_max[index] and fesom_y_hr[n] > y_min[index] and fesom_y_hr[n] < y_max[index]:
-                            x_index = nonzero(x_bins > fesom_x_hr[n])[0][0]-1
-                            y_index = nonzero(y_bins > fesom_y_hr[n])[0][0]-1
-                            fesom_u_hr[y_index, x_index] += u_circ_fesom_hr[n]
-                            fesom_v_hr[y_index, x_index] += v_circ_fesom_hr[n]
-                            fesom_num_pts_hr[y_index, x_index] += 1
+                    if fesom_x_hr[n] > x_min[index] and fesom_x_hr[n] < x_max[index] and fesom_y_hr[n] > y_min[index] and fesom_y_hr[n] < y_max[index]:
+                        x_index = nonzero(x_bins > fesom_x_hr[n])[0][0]-1
+                        y_index = nonzero(y_bins > fesom_y_hr[n])[0][0]-1
+                        fesom_u_hr[y_index, x_index] += u_circ_fesom_hr[n]
+                        fesom_v_hr[y_index, x_index] += v_circ_fesom_hr[n]
+                        fesom_num_pts_hr[y_index, x_index] += 1
                 fesom_u_hr = ma.masked_where(fesom_num_pts_hr==0, fesom_u_hr)
                 fesom_v_hr = ma.masked_where(fesom_num_pts_hr==0, fesom_v_hr)
                 flag = fesom_num_pts_hr > 0
@@ -552,16 +574,21 @@ def mip_regions_1var ():
             ax = fig.add_subplot(1,3,2, aspect='equal')
             # Start with land background
             contourf(x_reg_fesom, y_reg_fesom, land_square, 1, colors=(('0.6', '0.6', '0.6')))
-            # Add ice shelf elements
-            img = PatchCollection(patches_lr, cmap=colour_map)
+            # Add elements
+            if var in ['draft', 'melt', 'wct']:
+                # Ice shelf elements only
+                img = PatchCollection(patches_shelf_lr, cmap=colour_map)
+            else:
+                img = PatchCollection(patches_lr, cmap=colour_map)
             img.set_array(array(fesom_data_lr))
             img.set_edgecolor('face')
             img.set_clim(vmin=var_min, vmax=var_max)
             ax.add_collection(img)
-            # Mask out the open ocean in white
-            overlay = PatchCollection(mask_patches_lr, facecolor=(1,1,1))
-            overlay.set_edgecolor('face')
-            ax.add_collection(overlay)
+            if var in ['draft', 'melt', 'wct']:
+                # Mask out the open ocean in white
+                overlay = PatchCollection(patches_ocn_lr, facecolor=(1,1,1))
+                overlay.set_edgecolor('face')
+                ax.add_collection(overlay)
             if var == 'vel':
                 # Overlay vectors
                 quiver(x_centres, y_centres, fesom_u_lr, fesom_v_lr, scale=1.5, headwidth=6, headlength=7, color='black')
@@ -572,14 +599,19 @@ def mip_regions_1var ():
             # FESOM high-res
             ax = fig.add_subplot(1,3,3, aspect='equal')
             contourf(x_reg_fesom, y_reg_fesom, land_square, 1, colors=(('0.6', '0.6', '0.6')))
-            img = PatchCollection(patches_hr, cmap=colour_map)
+            if var in ['draft', 'melt', 'wct']:
+                # Ice shelf elements only
+                img = PatchCollection(patches_shelf_hr, cmap=colour_map)
+            else:
+                img = PatchCollection(patches_hr, cmap=colour_map)
             img.set_array(array(fesom_data_hr))
             img.set_edgecolor('face')
             img.set_clim(vmin=var_min, vmax=var_max)
             ax.add_collection(img)
-            overlay = PatchCollection(mask_patches_hr, facecolor=(1,1,1))
-            overlay.set_edgecolor('face')
-            ax.add_collection(overlay)
+            if var in ['draft', 'melt', 'wct']:
+                overlay = PatchCollection(patches_ocn_hr, facecolor=(1,1,1))
+                overlay.set_edgecolor('face')
+                ax.add_collection(overlay)
             if var == 'vel':
                 # Overlay vectors
                 quiver(x_centres, y_centres, fesom_u_hr, fesom_v_hr, scale=1.5, headwidth=6, headlength=7, color='black')
diff --git a/mip_scatterplot.py b/mip_scatterplot.py
index ae2eacc..6a135d1 100644
--- a/mip_scatterplot.py
+++ b/mip_scatterplot.py
@@ -1,7 +1,7 @@
 from numpy import *
 from matplotlib.pyplot import *
 
-def mip_scatterplot (roms_logfile, fesom_logfile_lr, fesom_logfile_hr):
+def mip_scatterplot (roms_logfile, roms_logfile_bs, fesom_logfile_lr, fesom_logfile_bs_lr, fesom_logfile_hr, fesom_logfile_bs_hr):
 
     # Year simulations start
     year_start = 1992
@@ -16,8 +16,15 @@ def mip_scatterplot (roms_logfile, fesom_logfile_lr, fesom_logfile_hr):
     obs_massloss = [1.4, 20.7, 135.4, 155.4, 51.8, 101.2, 97.5, 45.2, 144.9, 4.2, 18.2, 7.9, 90.6, 72.6, 27.2, 35.5, -2, 21.6, 6.3, 3.9, 26.8, 9.7, 47.7]
     obs_massloss_error = [14, 67, 40, 45, 19, 8, 7, 4, 14, 2, 3, 3, 8, 15, 10, 23, 3, 18, 2, 2, 14, 16, 34]
     num_shelves = len(obs_massloss)
-    # Order of indices for the ice shelves to be plotted on the x-axis (0-based)
-    order = [3, 21, 20, 19, 18, 17, 16, 15, 14, 13, 12, 11, 22, 10, 9, 8, 7, 6, 5, 4, 2, 1, 0]
+    # Some Bellingshausen ice shelves were split up later
+    names_bs = ['Wilkins Ice Shelf', 'Stange Ice Shelf', 'George VI Ice Shelf']
+    obs_massloss_bs = [18.4, 28, 89]
+    obs_massloss_error_bs = [17, 6, 17]
+    num_shelves_bs = len(obs_massloss_bs)
+    # Order of indices for the ice shelves to be plotted on the x-axis
+    # (0-based, assuming the Bellingshausen ice shelves have been tacked onto
+    # the end of the original arrays)
+    order = [3, 21, 20, 19, 18, 17, 16, 15, 14, 13, 12, 11, 22, 10, 9, 8, 7, 6, 5, 4, 24, 25, 23, 1, 0]
 
     # Read ROMS logfile
     roms_time = []
@@ -62,6 +69,35 @@ def mip_scatterplot (roms_logfile, fesom_logfile_lr, fesom_logfile_hr):
         t_end = nonzero(roms_time >= calc_end+1)[0][0]
     roms_massloss = mean(roms_massloss_ts[:,t_start:t_end], axis=1)
 
+    # Repeat for Bellingshausen
+    f = open(roms_logfile_bs, 'r')
+    f.readline()
+    # Skip the time values (should be the same)
+    for line in f:
+        try:
+            tmp = float(line)
+        except(ValueError):
+            # Reached the header for the next variable
+            break       
+    roms_massloss_bs_ts = empty([num_shelves_bs, len(roms_time)])
+    index = 0
+    while index < num_shelves_bs:
+        t = 0
+        for line in f:
+            try:
+                roms_massloss_bs_ts[index, t] = float(line)
+                t += 1
+            except(ValueError):
+                break
+        index +=1
+    f.close()
+    t_start = nonzero(roms_time >= calc_start)[0][0]
+    if calc_end == 2016:
+        t_end = size(roms_time)
+    else:
+        t_end = nonzero(roms_time >= calc_end+1)[0][0]
+    roms_massloss_bs = mean(roms_massloss_bs_ts[:,t_start:t_end], axis=1)
+
     # Read FESOM timeseries
     # Low-res
     f = open(fesom_logfile_lr, 'r')
@@ -93,6 +129,24 @@ def mip_scatterplot (roms_logfile, fesom_logfile_lr, fesom_logfile_hr):
     f.close()
     # Average between given years
     fesom_massloss_lr = mean(fesom_massloss_ts_lr[:,peryear*(calc_start-year_start):peryear*(calc_end+1-year_start)], axis=1)
+
+    # Repeat for Bellingshausen
+    f = open(fesom_logfile_bs_lr, 'r')
+    f.readline()
+    fesom_massloss_bs_ts_lr = empty([num_shelves_bs, num_time])
+    index = 0
+    while index < num_shelves_bs:
+        t = 0
+        for line in f:
+            try:
+                fesom_massloss_bs_ts_lr[index,t] = float(line)
+                t += 1
+            except(ValueError):
+                break
+        index += 1
+    f.close()
+    fesom_massloss_bs_lr = mean(fesom_massloss_bs_ts_lr[:,peryear*(calc_start-year_start):peryear*(calc_end+1-year_start)], axis=1)
+
     # Repeat for high-res
     f = open(fesom_logfile_hr, 'r')
     f.readline()
@@ -117,34 +171,65 @@ def mip_scatterplot (roms_logfile, fesom_logfile_lr, fesom_logfile_hr):
     f.close()
     fesom_massloss_hr = mean(fesom_massloss_ts_hr[:,peryear*(calc_start-year_start):peryear*(calc_end+1-year_start)], axis=1)
 
+    # High-res Bellingshausen
+    f = open(fesom_logfile_bs_hr, 'r')
+    f.readline()
+    fesom_massloss_bs_ts_hr = empty([num_shelves_bs, num_time])
+    index = 0
+    while index < num_shelves_bs:
+        t = 0
+        for line in f:
+            try:
+                fesom_massloss_bs_ts_hr[index,t] = float(line)
+                t += 1
+            except(ValueError):
+                break
+        index += 1
+    f.close()
+    fesom_massloss_bs_hr = mean(fesom_massloss_bs_ts_hr[:,peryear*(calc_start-year_start):peryear*(calc_end+1-year_start)], axis=1)
+
+    # Concatenate the Bellingshausen arrays onto the ends of the original arrays
+    names = names + names_bs
+    obs_massloss = obs_massloss + obs_massloss_bs
+    obs_massloss_error = obs_massloss_error + obs_massloss_error_bs
+    roms_massloss = concatenate((roms_massloss, roms_massloss_bs))
+    fesom_massloss_lr = concatenate((fesom_massloss_lr, fesom_massloss_bs_lr))
+    fesom_massloss_hr = concatenate((fesom_massloss_hr, fesom_massloss_bs_hr))
+    num_shelves_plot = len(order)
+
     # Figure out error values, in correct order for plotting
     roms_error = []
     fesom_error_lr = []
     fesom_error_hr = []
+    error_bars = []
     labels = []
-    for index in range(num_shelves):
-        obs_min = obs_massloss[order[index]] - obs_massloss_error[order[index]]
-        obs_max = obs_massloss[order[index]] + obs_massloss_error[order[index]]
-        if roms_massloss[order[index]] < obs_min:
-            roms_error.append(roms_massloss[order[index]] - obs_min)
-        elif roms_massloss[order[index]] > obs_max:
-            roms_error.append(roms_massloss[order[index]] - obs_max)
-        else:
-            roms_error.append(0)
-        if fesom_massloss_lr[order[index]] < obs_min:
-            fesom_error_lr.append(fesom_massloss_lr[order[index]] - obs_min)
-        elif fesom_massloss_lr[order[index]] > obs_max:
-            fesom_error_lr.append(fesom_massloss_lr[order[index]] - obs_max)
-        else:
-            fesom_error_lr.append(0)
-        if fesom_massloss_hr[order[index]] < obs_min:
-            fesom_error_hr.append(fesom_massloss_hr[order[index]] - obs_min)
-        elif fesom_massloss_hr[order[index]] > obs_max:
-            fesom_error_hr.append(fesom_massloss_hr[order[index]] - obs_max)
-        else:
-            fesom_error_hr.append(0)
-        labels.append(names[order[index]])
-
+    for index in order:
+        #obs_min = obs_massloss[index] - obs_massloss_error[index]
+        #obs_max = obs_massloss[index] + obs_massloss_error[index]
+        #if roms_massloss[index] < obs_min:
+        #    roms_error.append(roms_massloss[index] - obs_min)
+        #elif roms_massloss[index] > obs_max:
+        #    roms_error.append(roms_massloss[index] - obs_max)
+        #else:
+        #    roms_error.append(0)
+        #if fesom_massloss_lr[index] < obs_min:
+        #    fesom_error_lr.append(fesom_massloss_lr[index] - obs_min)
+        #elif fesom_massloss_lr[index] > obs_max:
+        #    fesom_error_lr.append(fesom_massloss_lr[index] - obs_max)
+        #else:
+        #    fesom_error_lr.append(0)
+        #if fesom_massloss_hr[index] < obs_min:
+        #    fesom_error_hr.append(fesom_massloss_hr[index] - obs_min)
+        #elif fesom_massloss_hr[index] > obs_max:
+        #    fesom_error_hr.append(fesom_massloss_hr[index] - obs_max)
+        #else:
+        #    fesom_error_hr.append(0)
+        roms_error.append(roms_massloss[index] - obs_massloss[index])
+        fesom_error_lr.append(fesom_massloss_lr[index] - obs_massloss[index])
+        fesom_error_hr.append(fesom_massloss_hr[index] - obs_massloss[index])
+        error_bars.append(obs_massloss_error[index])
+        labels.append(names[index])
+    
     # Plot
     fig = figure(figsize=(10,7))
     gs = GridSpec(1,1)
@@ -154,28 +239,31 @@ def mip_scatterplot (roms_logfile, fesom_logfile_lr, fesom_logfile_hr):
     axvspan(0.5, 6.5, facecolor='b', alpha=0.1)
     axvspan(7.5, 11.5, facecolor='b', alpha=0.1)
     axvspan(14.5, 18.5, facecolor='b', alpha=0.1)
-    axvspan(20.5, 23, facecolor='b', alpha=0.1)
+    axvspan(22.5, 25, facecolor='b', alpha=0.1)
     # Region labels
-    text(0, 30, 'FR', fontsize=14, ha='center')
-    text(3.5, 30, 'EWed', fontsize=14, ha='center')
-    text(7, -20, 'Am', fontsize=14, ha='center')
-    text(9.5, 30, 'Aus', fontsize=14, ha='center')
-    text(13, 30, 'RS', fontsize=14, ha='center')
-    text(16.5, 30, 'AS', fontsize=14, ha='center')
-    text(19.5, 30, 'BS', fontsize=14, ha='center')
-    text(21.5, 30, 'Lr', fontsize=14, ha='center')
-    plot(range(-1,num_shelves+1), zeros(num_shelves+2), color='black', linewidth=3)
-    plot(range(num_shelves), roms_error, 'o', color=(0.08, 0.4, 0.79), ms=10, label='MetROMS')
-    plot(range(num_shelves), fesom_error_lr, 'o', color=(0.73, 0.06, 0.69), ms=10, label='FESOM (low-res)')
-    plot(range(num_shelves), fesom_error_hr, 'o', color=(0.06, 0.73, 0.1), ms=10, label='FESOM (high-res)')
+    text(0, 80, 'FR', fontsize=14, ha='center')
+    text(3.5, 80, 'EWed', fontsize=14, ha='center')
+    text(7, 80, 'Am', fontsize=14, ha='center')
+    text(9.5, 80, 'Aus', fontsize=14, ha='center')
+    text(13, 80, 'RS', fontsize=14, ha='center')
+    text(16.5, 80, 'AS', fontsize=14, ha='center')
+    text(20.5, 80, 'BS', fontsize=14, ha='center')
+    text(23.5, 80, 'Lr', fontsize=14, ha='center')
+    # Black line at zero
+    plot(range(-1,num_shelves_plot+1), zeros(num_shelves_plot+2), color='black', linewidth=3)
+    plot(range(num_shelves_plot), roms_error, 'o', color=(0.08, 0.4, 0.79), ms=10, label='MetROMS')
+    plot(range(num_shelves_plot), fesom_error_lr, 'o', color=(0.06, 0.73, 0.1), ms=10, label='FESOM (low-res)')
+    plot(range(num_shelves_plot), fesom_error_hr, 'o', color=(0.73, 0.06, 0.69), ms=10, label='FESOM (high-res)')
     grid(True)    
-    xlim([-0.5, num_shelves-0.5])
-    xticks(range(num_shelves), labels, rotation=90)
+    xlim([-0.5, num_shelves_plot-0.5])
+    xticks(range(num_shelves_plot), labels, rotation=90)
     ylabel('Gt/y', fontsize=14)
     title('Bias in Ice Shelf Basal Mass Loss', fontsize=18)
     # Make legend
-    legend(numpoints=1,loc='lower left')
+    legend(numpoints=1,bbox_to_anchor=(0.75,-0.3)) #loc='lower left')
     setp(gca().get_legend().get_texts(), fontsize='13')
+    # Error bars on top
+    errorbar(range(num_shelves_plot), zeros(num_shelves_plot), yerr=error_bars, color='black', capthick=1)
     fig.show()
     fig.savefig('scatterplot.png')
 
@@ -184,7 +272,9 @@ def mip_scatterplot (roms_logfile, fesom_logfile_lr, fesom_logfile_hr):
 if __name__ == "__main__":
 
     roms_logfile = raw_input("Path to ROMS logfile from timeseries_massloss.py: ")
+    roms_logfile_bs = raw_input("Path to ROMS logfile from timeseries_massloss_bellingshausen.py: ")
     fesom_logfile_lr = raw_input("Path to FESOM low-res logfile from timeseries_massloss.py: ")
+    fesom_logfile_bs_lr = raw_input("Path to FESOM low-res logfile from timeseries_massloss_bellingshausen.py: ")
     fesom_logfile_hr = raw_input("Path to FESOM high-res logfile from timeseries_massloss.py: ")
-    mip_scatterplot(roms_logfile, fesom_logfile_lr, fesom_logfile_hr)
-
+    fesom_logfile_bs_hr = raw_input("Path to FESOM high-res logfile from timeseries_massloss_bellingshausen.py: ")
+    mip_scatterplot(roms_logfile, roms_logfile_bs, fesom_logfile_lr, fesom_logfile_bs_lr, fesom_logfile_hr, fesom_logfile_bs_hr)
diff --git a/mip_seasonal_cycle.py b/mip_seasonal_cycle.py
index 5e62cdc..31233dd 100644
--- a/mip_seasonal_cycle.py
+++ b/mip_seasonal_cycle.py
@@ -1,8 +1,8 @@
 from numpy import *
 
 # Calculate the average amplitude of the seasonal cycle in total basal mass
-# loss over 2002-2016, in MetROMS, low-res FESOM, and high-res FESOM. Print
-# results to the screen.
+# loss over 2002-2016, in MetROMS, low-res FESOM, and high-res FESOM, for the
+# entire continent as well as the Amery. Print results to the screen.
 def mip_seasonal_cycle (roms_logfile, fesom_logfile_lr, fesom_logfile_hr):
 
     # Year simulations start
@@ -12,6 +12,11 @@ def mip_seasonal_cycle (roms_logfile, fesom_logfile_lr, fesom_logfile_hr):
     calc_end = 2016
     # Number of output steps per year in FESOM
     peryear = 365/5
+    # Name of each ice shelf
+    names = ['Total Mass Loss', 'Larsen D Ice Shelf', 'Larsen C Ice Shelf', 'Wilkins & George VI & Stange Ice Shelves', 'Ronne-Filchner Ice Shelf', 'Abbot Ice Shelf', 'Pine Island Glacier Ice Shelf', 'Thwaites Ice Shelf', 'Dotson Ice Shelf', 'Getz Ice Shelf', 'Nickerson Ice Shelf', 'Sulzberger Ice Shelf', 'Mertz Ice Shelf', 'Totten & Moscow University Ice Shelves', 'Shackleton Ice Shelf', 'West Ice Shelf', 'Amery Ice Shelf', 'Prince Harald Ice Shelf', 'Baudouin & Borchgrevink Ice Shelves', 'Lazarev Ice Shelf', 'Nivl Ice Shelf', 'Fimbul & Jelbart & Ekstrom Ice Shelves', 'Brunt & Riiser-Larsen Ice Shelves', 'Ross Ice Shelf']
+    num_shelves = len(names)
+    i_total = 0
+    i_amery = 16
 
     # Read ROMS logfile
     roms_time = []
@@ -24,19 +29,30 @@ def mip_seasonal_cycle (roms_logfile, fesom_logfile_lr, fesom_logfile_hr):
         except(ValueError):
             # Reached the header for the next variable
             break
-    # Read total mass loss
-    roms_massloss = []
-    for line in f:
-        try:
-            roms_massloss.append(float(line))
-        except(ValueError):
-            # Reached the header for the first individual ice shelf
-            break
+    # Set up array for mass loss values at each ice shelf
+    roms_massloss = empty([num_shelves, len(roms_time)])
+    index = 0
+    # Loop over ice shelves
+    while index < num_shelves:
+        t = 0
+        for line in f:
+            try:
+                roms_massloss[index, t] = float(line)
+                t += 1
+            except(ValueError):
+                # Reached the header for the next ice shelf
+                break
+        index +=1
     f.close()
     # Add start year to ROMS time array
     roms_time = array(roms_time) + year_start
     # Calculate amplitude for each year
-    roms_amplitude = []
+    #roms_amplitude_total = []
+    roms_min_total = []
+    roms_max_total = []
+    #roms_amplitude_amery = []
+    roms_min_amery = []
+    roms_max_amery = []
     for year in range(calc_start, calc_end):
         # Find the first index after the beginning of this year
         t_start = nonzero(roms_time >= year)[0][0]
@@ -48,15 +64,29 @@ def mip_seasonal_cycle (roms_logfile, fesom_logfile_lr, fesom_logfile_hr):
             t_end = len(roms_time)
         else:
             t_end = tmp2[0]
+        # Total Antarctica
         # Get min and max between these bounds
-        curr_min = amin(roms_massloss[t_start:t_end])
-        curr_max = amax(roms_massloss[t_start:t_end])
+        curr_min = amin(roms_massloss[i_total,t_start:t_end])
+        curr_max = amax(roms_massloss[i_total,t_start:t_end])
+        roms_min_total.append(curr_min)
+        roms_max_total.append(curr_max)
         # Save amplitude
-        roms_amplitude.append(curr_max-curr_min)
+        #roms_amplitude_total.append(curr_max-curr_min)
+        # Amery
+        curr_min = amin(roms_massloss[i_amery,t_start:t_end])
+        curr_max = amax(roms_massloss[i_amery,t_start:t_end])
+        #roms_amplitude_amery.append(curr_max-curr_min)
+        roms_min_amery.append(curr_min)
+        roms_max_amery.append(curr_max)
     # Calculate average amplitude
-    roms_val = mean(array(roms_amplitude))
-    print 'Average amplitude of seasonal cycle in total basal mass loss: '
-    print 'ROMS: ' + str(roms_val)
+    #roms_val_total = mean(array(roms_amplitude_total))
+    #roms_val_amery = mean(array(roms_amplitude_amery))
+    #print 'Average amplitude of seasonal cycle in total basal mass loss: '
+    #print 'ROMS (Total Antarctica): ' + str(roms_val_total)
+    #print 'ROMS (Amery): ' + str(roms_val_amery)
+    print 'Average min/max in total basal mass loss: '
+    print 'ROMS (Total Antarctica): ' + str(mean(array(roms_min_total))) + ' Gt/y min, ' + str(mean(array(roms_max_total))) + ' Gt/y max'
+    print 'ROMS (Amery): ' + str(mean(array(roms_min_amery))) + ' Gt/y min, ' + str(mean(array(roms_max_amery))) + ' Gt/y max'
 
     # Read FESOM logfiles
     # Low-res
@@ -64,47 +94,113 @@ def mip_seasonal_cycle (roms_logfile, fesom_logfile_lr, fesom_logfile_hr):
     # Skip the first line (header)
     f.readline()
     # Read total mass loss
-    fesom_massloss_lr = []
+    tmp = []
+    num_time = 0
     for line in f:
         try:
-            fesom_massloss_lr.append(float(line))
+            tmp.append(float(line))
+            num_time += 1
         except(ValueError):
             # Reached the header for the first individual ice shelf
             break
+    # Set up array for mass loss values at each ice shelf
+    fesom_massloss_lr = empty([num_shelves, num_time])
+    # Save the total values in the first index
+    fesom_massloss_lr[0,:] = array(tmp)
+    # Loop over ice shelves
+    index = 1
+    while index < num_shelves:
+        t = 0
+        for line in f:
+            try:
+                fesom_massloss_lr[index,t] = float(line)
+                t += 1
+            except(ValueError):
+                # Reached the header for the next ice shelf
+                break
+        index += 1
     f.close()
     # Calculate amplitude for each year
-    fesom_amplitude_lr = []
+    #fesom_amplitude_total_lr = []
+    fesom_min_total_lr = []
+    fesom_max_total_lr = []
+    #fesom_amplitude_amery_lr = []
+    fesom_min_amery_lr = []
+    fesom_max_amery_lr = []
     for year in range(calc_start, calc_end):
         t_start = peryear*(year-year_start)
         t_end = peryear*(year+1-year_start)
+        # Total Antarctica
         # Get min and max between these bounds
-        curr_min = amin(fesom_massloss_lr[t_start:t_end])
-        curr_max = amax(fesom_massloss_lr[t_start:t_end])
+        curr_min = amin(fesom_massloss_lr[i_total,t_start:t_end])
+        curr_max = amax(fesom_massloss_lr[i_total,t_start:t_end])
         # Save amplitude
-        fesom_amplitude_lr.append(curr_max-curr_min)
+        #fesom_amplitude_total_lr.append(curr_max-curr_min)
+        fesom_min_total_lr.append(curr_min)
+        fesom_max_total_lr.append(curr_max)
+        # Amery
+        curr_min = amin(fesom_massloss_lr[i_amery,t_start:t_end])
+        curr_max = amax(fesom_massloss_lr[i_amery,t_start:t_end])
+        #fesom_amplitude_amery_lr.append(curr_max-curr_min)
+        fesom_min_amery_lr.append(curr_min)
+        fesom_max_amery_lr.append(curr_max)
     # Calculate average amplitude
-    fesom_val_lr = mean(array(fesom_amplitude_lr))
-    print 'FESOM low-res: ' + str(fesom_val_lr)
+    #fesom_val_total_lr = mean(array(fesom_amplitude_total_lr))
+    #fesom_val_amery_lr = mean(array(fesom_amplitude_amery_lr))
+    #print 'FESOM low-res (Total Antarctica): ' + str(fesom_val_total_lr)
+    #print 'FESOM low-res (Amery): ' + str(fesom_val_amery_lr)
+    print 'FESOM low-res (Total Antarctica): ' + str(mean(array(fesom_min_total_lr))) + ' Gt/y min, ' + str(mean(array(fesom_max_total_lr))) + ' Gt/y max'
+    print 'FESOM low-res (Amery): ' + str(mean(array(fesom_min_amery_lr))) + ' Gt/y min, ' + str(mean(array(fesom_max_amery_lr))) + ' Gt/y max'
 
     # High-res
     f = open(fesom_logfile_hr, 'r')
     f.readline()
-    fesom_massloss_hr = []
+    tmp = []
+    num_time = 0
     for line in f:
         try:
-            fesom_massloss_hr.append(float(line))
+            tmp.append(float(line))
+            num_time += 1
         except(ValueError):
             break
+    fesom_massloss_hr = empty([num_shelves, num_time])
+    fesom_massloss_hr[0,:] = array(tmp)
+    index = 1
+    while index < num_shelves:
+        t = 0
+        for line in f:
+            try:
+                fesom_massloss_hr[index,t] = float(line)
+                t += 1
+            except(ValueError):
+                break
+        index += 1
     f.close()
-    fesom_amplitude_hr = []
+    #fesom_amplitude_total_hr = []
+    fesom_min_total_hr = []
+    fesom_max_total_hr = []
+    #fesom_amplitude_amery_hr = []
+    fesom_min_amery_hr = []
+    fesom_max_amery_hr = []
     for year in range(calc_start, calc_end):
         t_start = peryear*(year-year_start)
         t_end = peryear*(year+1-year_start)
-        curr_min = amin(fesom_massloss_hr[t_start:t_end])
-        curr_max = amax(fesom_massloss_hr[t_start:t_end])
-        fesom_amplitude_hr.append(curr_max-curr_min)
-    fesom_val_hr = mean(array(fesom_amplitude_hr))
-    print 'FESOM high-res: ' + str(fesom_val_hr)
+        curr_min = amin(fesom_massloss_hr[i_total,t_start:t_end])
+        curr_max = amax(fesom_massloss_hr[i_total,t_start:t_end])
+        fesom_min_total_hr.append(curr_min)
+        fesom_max_total_hr.append(curr_max)
+        #fesom_amplitude_total_hr.append(curr_max-curr_min)
+        curr_min = amin(fesom_massloss_hr[i_amery,t_start:t_end])
+        curr_max = amax(fesom_massloss_hr[i_amery,t_start:t_end])
+        #fesom_amplitude_amery_hr.append(curr_max-curr_min)
+        fesom_min_amery_hr.append(curr_min)
+        fesom_max_amery_hr.append(curr_max)
+    #fesom_val_total_hr = mean(array(fesom_amplitude_total_hr))
+    #fesom_val_amery_hr = mean(array(fesom_amplitude_amery_hr))
+    #print 'FESOM high-res (Total Antarctica): ' + str(fesom_val_total_hr)
+    #print 'FESOM high-res (Amery): ' + str(fesom_val_amery_hr)
+    print 'FESOM high-res (Total Antarctica): ' + str(mean(array(fesom_min_total_hr))) + ' Gt/y min, ' + str(mean(array(fesom_max_total_hr))) + ' Gt/y max'
+    print 'FESOM high-res (Amery): ' + str(mean(array(fesom_min_amery_hr))) + ' Gt/y min, ' + str(mean(array(fesom_max_amery_hr))) + ' Gt/y max'
         
 
 # Command-line interface
diff --git a/mip_ts_distribution.py b/mip_ts_distribution.py
index 72ae325..4d732a0 100644
--- a/mip_ts_distribution.py
+++ b/mip_ts_distribution.py
@@ -16,10 +16,12 @@
 # roms_grid = path to ROMS grid file
 # roms_file = path to time-averaged ROMS file containing temperature and
 #             salinity (I used 2002-2016 average)
-# fesom_mesh_path = path to FESOM mesh directory (I used high-res)
-# fesom_file = path to time-averaged FESOM file containing temperature and
-#              salinity, over the same period as roms_file
-def mip_ts_distribution (roms_grid, roms_file, fesom_mesh_path, fesom_file):
+# fesom_mesh_path_lr, fesom_mesh_path_hr = path to FESOM mesh directories for
+#                     low-res and high-res meshes
+# fesom_file_lr, fesom_file_hr = paths to time-averaged FESOM files containing
+#                temperature and salinity for the low-res and high-res
+#                simulations respectively, over the same period as roms_file
+def mip_ts_distribution (roms_grid, roms_file, fesom_mesh_path_lr, fesom_file_lr, fesom_mesh_path_hr, fesom_file_hr):
 
     # Northern boundary of water masses to consider
     nbdry = -65
@@ -32,7 +34,7 @@ def mip_ts_distribution (roms_grid, roms_file, fesom_mesh_path, fesom_file):
     max_temp = 3.8
     # Bounds to actually plot
     min_salt_plot = 33.25
-    max_salt_plot = 35.0
+    max_salt_plot = 35.1
     min_temp_plot = -3
     max_temp_plot = 3.8
     # FESOM grid generation parameters
@@ -55,10 +57,12 @@ def mip_ts_distribution (roms_grid, roms_file, fesom_mesh_path, fesom_file):
     # Set up 2D arrays of temperature bins x salinity bins to hold average
     # depth of water masses, weighted by volume
     ts_vals_roms = zeros([size(temp_centres), size(salt_centres)])
-    ts_vals_fesom = zeros([size(temp_centres), size(salt_centres)])
+    ts_vals_fesom_lr = zeros([size(temp_centres), size(salt_centres)])
+    ts_vals_fesom_hr = zeros([size(temp_centres), size(salt_centres)])
     # Also arrays to integrate volume
     volume_roms = zeros([size(temp_centres), size(salt_centres)])
-    volume_fesom = zeros([size(temp_centres), size(salt_centres)])
+    volume_fesom_lr = zeros([size(temp_centres), size(salt_centres)])
+    volume_fesom_hr = zeros([size(temp_centres), size(salt_centres)])
     # Calculate surface freezing point as a function of salinity as seen by
     # each sea ice model
     freezing_pt_roms = salt_centres/(-18.48 + 18.48/1e3*salt_centres)
@@ -112,16 +116,16 @@ def mip_ts_distribution (roms_grid, roms_file, fesom_mesh_path, fesom_file):
     # Convert depths from integrals to volume-averages
     ts_vals_roms /= volume_roms
 
-    print 'Processing FESOM'
+    print 'Processing low-res FESOM'
     # Make FESOM grid elements
-    elements = fesom_grid(fesom_mesh_path, circumpolar, cross_180)
+    elements_lr = fesom_grid(fesom_mesh_path_lr, circumpolar, cross_180)
     # Read temperature and salinity at each 3D node
-    id = Dataset(fesom_file, 'r')
-    fesom_temp = id.variables['temp'][0,:]
-    fesom_salt = id.variables['salt'][0,:]
+    id = Dataset(fesom_file_lr, 'r')
+    fesom_temp_lr = id.variables['temp'][0,:]
+    fesom_salt_lr = id.variables['salt'][0,:]
     id.close()
     # Loop over elements
-    for elm in elements:
+    for elm in elements_lr:
         # See if we're in the region of interest
         if all(elm.lat < nbdry):
             # Get area of 2D triangle
@@ -140,11 +144,11 @@ def mip_ts_distribution (roms_grid, roms_file, fesom_mesh_path, fesom_file):
                 dz = []
                 for i in range(3):
                     # Average temperature over 6 nodes
-                    temp_vals.append(fesom_temp[nodes[i].id])
-                    temp_vals.append(fesom_temp[nodes[i].below.id])
+                    temp_vals.append(fesom_temp_lr[nodes[i].id])
+                    temp_vals.append(fesom_temp_lr[nodes[i].below.id])
                     # Average salinity over 6 nodes
-                    salt_vals.append(fesom_salt[nodes[i].id])
-                    salt_vals.append(fesom_salt[nodes[i].below.id])
+                    salt_vals.append(fesom_salt_lr[nodes[i].id])
+                    salt_vals.append(fesom_salt_lr[nodes[i].below.id])
                     # Average depth over 6 nodes
                     depth_vals.append(nodes[i].depth)
                     depth_vals.append(nodes[i].below.depth)
@@ -161,26 +165,64 @@ def mip_ts_distribution (roms_grid, roms_file, fesom_mesh_path, fesom_file):
                 temp_index = nonzero(temp_bins > temp_elm)[0][0] - 1
                 salt_index = nonzero(salt_bins > salt_elm)[0][0] - 1
                 # Integrate depth*volume in this bin
-                ts_vals_fesom[temp_index, salt_index] += depth_elm*volume
-                volume_fesom[temp_index, salt_index] += volume
+                ts_vals_fesom_lr[temp_index, salt_index] += depth_elm*volume
+                volume_fesom_lr[temp_index, salt_index] += volume
     # Mask bins with zero volume
-    ts_vals_fesom = ma.masked_where(volume_fesom==0, ts_vals_fesom)
-    volume_fesom = ma.masked_where(volume_fesom==0, volume_fesom)
+    ts_vals_fesom_lr = ma.masked_where(volume_fesom_lr==0, ts_vals_fesom_lr)
+    volume_fesom_lr = ma.masked_where(volume_fesom_lr==0, volume_fesom_lr)
     # Convert depths from integrals to volume-averages
-    ts_vals_fesom /= volume_fesom
+    ts_vals_fesom_lr /= volume_fesom_lr
+
+    print 'Processing high-res FESOM'
+    elements_hr = fesom_grid(fesom_mesh_path_hr, circumpolar, cross_180)
+    id = Dataset(fesom_file_hr, 'r')
+    fesom_temp_hr = id.variables['temp'][0,:]
+    fesom_salt_hr = id.variables['salt'][0,:]
+    id.close()
+    for elm in elements_hr:
+        if all(elm.lat < nbdry):
+            area = elm.area()
+            nodes = [elm.nodes[0], elm.nodes[1], elm.nodes[2]]
+            while True:
+                if nodes[0].below is None or nodes[1].below is None or nodes[2].below is None:
+                    break
+                temp_vals = []
+                salt_vals = []
+                depth_vals = []
+                dz = []
+                for i in range(3):
+                    temp_vals.append(fesom_temp_hr[nodes[i].id])
+                    temp_vals.append(fesom_temp_hr[nodes[i].below.id])
+                    salt_vals.append(fesom_salt_hr[nodes[i].id])
+                    salt_vals.append(fesom_salt_hr[nodes[i].below.id])
+                    depth_vals.append(nodes[i].depth)
+                    depth_vals.append(nodes[i].below.depth)
+                    dz.append(abs(nodes[i].depth - nodes[i].below.depth))
+                    nodes[i] = nodes[i].below
+                temp_elm = mean(array(temp_vals))
+                salt_elm = mean(array(salt_vals))
+                depth_elm = mean(array(depth_vals))
+                volume = area*mean(array(dz))
+                temp_index = nonzero(temp_bins > temp_elm)[0][0] - 1
+                salt_index = nonzero(salt_bins > salt_elm)[0][0] - 1
+                ts_vals_fesom_hr[temp_index, salt_index] += depth_elm*volume
+                volume_fesom_hr[temp_index, salt_index] += volume
+    ts_vals_fesom_hr = ma.masked_where(volume_fesom_hr==0, ts_vals_fesom_hr)
+    volume_fesom_hr = ma.masked_where(volume_fesom_hr==0, volume_fesom_hr)
+    ts_vals_fesom_hr /= volume_fesom_hr
 
     # Find the maximum depth for plotting
-    max_depth = max(amax(ts_vals_roms), amax(ts_vals_fesom))
+    max_depth = amax(array([amax(ts_vals_roms), amax(ts_vals_fesom_lr), amax(ts_vals_fesom_hr)]))
     # Make a nonlinear scale
     bounds = linspace(0, max_depth**(1.0/2.5), num=100)**2.5
     norm = BoundaryNorm(boundaries=bounds, ncolors=256)
     # Set labels for density contours
-    manual_locations = [(33.4, 3.0), (33.65, 3.0), (33.9, 3.0), (34.2, 3.0), (34.45, 3.5), (34.65, 3.25), (34.9, 3.0), (34.8, 0)]
+    manual_locations = [(33.4, 3.0), (33.65, 3.0), (33.9, 3.0), (34.2, 3.0), (34.45, 3.5), (34.65, 3.25), (34.9, 3.0), (35, 1.5)]
 
     print "Plotting"
-    fig = figure(figsize=(20,12))
+    fig = figure(figsize=(20,9))
     # ROMS
-    ax = fig.add_subplot(1, 2, 1)
+    ax = fig.add_subplot(1, 3, 1)
     pcolor(salt_centres, temp_centres, ts_vals_roms, norm=norm, vmin=0, vmax=max_depth, cmap='jet')
     # Add surface freezing point line
     plot(salt_centres, freezing_pt_roms, color='black', linestyle='dashed')
@@ -191,12 +233,24 @@ def mip_ts_distribution (roms_grid, roms_file, fesom_mesh_path, fesom_file):
     ylim([min_temp_plot, max_temp_plot])
     ax.tick_params(axis='x', labelsize=16)
     ax.tick_params(axis='y', labelsize=16)
-    xlabel('Salinity (psu)', fontsize=22)
-    ylabel(r'Temperature ($^{\circ}$C)', fontsize=22)
-    title('MetROMS', fontsize=28)
-    # FESOM
-    ax = fig.add_subplot(1, 2, 2)
-    img = pcolor(salt_centres, temp_centres, ts_vals_fesom, norm=norm, vmin=0, vmax=max_depth, cmap='jet')
+    xlabel('Salinity (psu)', fontsize=20)
+    ylabel(r'Temperature ($^{\circ}$C)', fontsize=20)
+    title('MetROMS', fontsize=24)
+    # FESOM low-res
+    ax = fig.add_subplot(1, 3, 2)
+    img = pcolor(salt_centres, temp_centres, ts_vals_fesom_lr, norm=norm, vmin=0, vmax=max_depth, cmap='jet')
+    plot(salt_centres, freezing_pt_fesom, color='black', linestyle='dashed')
+    cs = contour(salt_centres, temp_centres, density, density_lev, colors=(0.6,0.6,0.6), linestyles='dotted')
+    clabel(cs, inline=1, fontsize=14, color=(0.6,0.6,0.6), fmt='%1.1f', manual=manual_locations)
+    xlim([min_salt_plot, max_salt_plot])
+    ylim([min_temp_plot, max_temp_plot])
+    ax.tick_params(axis='x', labelsize=16)
+    ax.tick_params(axis='y', labelsize=16)
+    xlabel('Salinity (psu)', fontsize=20)
+    title('FESOM (low-res)', fontsize=24)
+    # FESOM high-res
+    ax = fig.add_subplot(1, 3, 3)
+    img = pcolor(salt_centres, temp_centres, ts_vals_fesom_hr, norm=norm, vmin=0, vmax=max_depth, cmap='jet')
     plot(salt_centres, freezing_pt_fesom, color='black', linestyle='dashed')
     cs = contour(salt_centres, temp_centres, density, density_lev, colors=(0.6,0.6,0.6), linestyles='dotted')
     clabel(cs, inline=1, fontsize=14, color=(0.6,0.6,0.6), fmt='%1.1f', manual=manual_locations)
@@ -204,8 +258,8 @@ def mip_ts_distribution (roms_grid, roms_file, fesom_mesh_path, fesom_file):
     ylim([min_temp_plot, max_temp_plot])
     ax.tick_params(axis='x', labelsize=16)
     ax.tick_params(axis='y', labelsize=16)
-    xlabel('Salinity (psu)', fontsize=22)
-    title('FESOM (high-res)', fontsize=28)
+    xlabel('Salinity (psu)', fontsize=20)
+    title('FESOM (high-res)', fontsize=24)
     # Add a colourbar on the right
     cbaxes = fig.add_axes([0.93, 0.2, 0.02, 0.6])
     cbar = colorbar(img, cax=cbaxes, ticks=[0,50,100,200,500,1000,2000,4000])
@@ -213,7 +267,7 @@ def mip_ts_distribution (roms_grid, roms_file, fesom_mesh_path, fesom_file):
     # Add the main title
     suptitle('Water masses south of 65$^{\circ}$S: depth (m), 2002-2016 average', fontsize=30)
     subplots_adjust(wspace=0.1)
-    #fig.show()
+    fig.show()
     fig.savefig('ts_distribution_orig.png')
 
 
@@ -222,9 +276,11 @@ def mip_ts_distribution (roms_grid, roms_file, fesom_mesh_path, fesom_file):
 
     roms_grid = raw_input("Path to ROMS grid file: ")
     roms_file = raw_input("Path to time-averaged ROMS file containing temperature and salinity: ")
-    fesom_mesh_path = raw_input("Path to FESOM mesh directory: ")
-    fesom_file = raw_input("Path to time-averaged FESOM file containing temperature and salinity: ")
-    mip_ts_distribution(roms_grid, roms_file, fesom_mesh_path, fesom_file)
+    fesom_mesh_path_lr = raw_input("Path to FESOM low-res mesh directory: ")
+    fesom_file_lr = raw_input("Path to time-averaged FESOM low-res file containing temperature and salinity: ")
+    fesom_mesh_path_hr = raw_input("Path to FESOM high-res mesh directory: ")
+    fesom_file_hr = raw_input("Path to time-averaged FESOM high-res file containing temperature and salinity: ")
+    mip_ts_distribution(roms_grid, roms_file, fesom_mesh_path_lr, fesom_file_lr, fesom_mesh_path_hr, fesom_file_hr)
     
             
       
diff --git a/mip_watermass_barchart.py b/mip_watermass_barchart.py
new file mode 100644
index 0000000..9de91bd
--- /dev/null
+++ b/mip_watermass_barchart.py
@@ -0,0 +1,321 @@
+from netCDF4 import Dataset
+from numpy import *
+from matplotlib.pyplot import *
+from cartesian_grid_3d import *
+import sys
+sys.path.insert(0, '/short/y99/kaa561/fesomtools')
+from fesom_grid import *
+
+def mip_watermass_barchart (roms_grid, roms_file, fesom_mesh_lr, fesom_mesh_hr, fesom_file_lr, fesom_file_hr):
+
+    # Sectors to consider
+    sector_names = ['a) Filchner-Ronne Ice Shelf', 'b) Eastern Weddell Region', 'c) Amery Ice Shelf', 'd) Australian Sector', 'e) Ross Sea', 'f) Amundsen Sea', 'g) Bellingshausen Sea', 'h) Larsen Ice Shelves', 'i) All Ice Shelves']
+    num_sectors = len(sector_names)
+    # Water masses to consider
+    wm_names = ['ISW', 'MCDW', 'HSSW', 'LSSW', 'AASW']
+    num_watermasses = len(wm_names)
+    wm_colours = [(0.73, 0.6, 1), (1, 0.4, 0.4), (0.52, 0.88, 0.52), (0.6, 0.8, 1), (1, 1, 0)]
+    # ROMS vertical grid parameters
+    theta_s = 7.0
+    theta_b = 2.0
+    hc = 250
+    N = 31
+    # FESOM mesh parameters
+    circumpolar = True
+    cross_180 = False
+
+    print 'Processing MetROMS'
+    # Read ROMS grid variables we need
+    id = Dataset(roms_grid, 'r')
+    roms_lon = id.variables['lon_rho'][:,:]
+    roms_lat = id.variables['lat_rho'][:,:]
+    roms_h = id.variables['h'][:,:]
+    roms_zice = id.variables['zice'][:,:]
+    id.close()
+    num_lat = size(roms_lat, 0)
+    num_lon = size(roms_lon, 1)
+    # Get integrands on 3D grid
+    roms_dx, roms_dy, roms_dz, roms_z = cartesian_grid_3d(roms_lon, roms_lat, roms_h, roms_zice, theta_s, theta_b, hc, N)
+    # Get volume integrand
+    dV = roms_dx*roms_dy*roms_dz
+    # Read ROMS output
+    id = Dataset(roms_file, 'r')
+    roms_temp = id.variables['temp'][0,:,:,:]
+    roms_salt = id.variables['salt'][0,:,:,:]
+    id.close()
+    # Initialise volume of each water mass in each sector
+    roms_vol_watermass = zeros([num_watermasses, num_sectors])
+    # Calculate water mass breakdown
+    for j in range(num_lat):
+        for i in range(num_lon):
+            # Select ice shelf points
+            if roms_zice[j,i] < 0:
+                # Figure out which sector this point falls into
+                lon = roms_lon[j,i]
+                if lon > 180:
+                    lon -= 360
+                lat = roms_lat[j,i]
+                if lon >= -85 and lon < -30 and lat < -74:
+                    # Filchner-Ronne
+                    sector = 0
+                elif lon >= -30 and lon < 65:
+                    # Eastern Weddell region
+                    sector = 1
+                elif lon >= 65 and lon < 76:
+                    # Amery
+                    sector = 2
+                elif lon >= 76 and lon < 165 and lat >= -74:
+                    # Australian sector
+                    sector = 3
+                elif (lon >= 155 and lon < 165 and lat < -74) or (lon >= 165) or (lon < -140):
+                    # Ross Sea
+                    sector = 4
+                elif (lon >= -140 and lon < -105) or (lon >= -105 and lon < -98 and lat < -73.1):
+                    # Amundsen Sea
+                    sector = 5
+                elif (lon >= -104 and lon < -98 and lat >= -73.1) or (lon >= -98 and lon < -66 and lat >= -75):
+                    # Bellingshausen Sea
+                    sector = 6
+                elif lon >= -66 and lon < -59 and lat >= -74:
+                    # Larsen Ice Shelves
+                    sector = 7
+                else:
+                    print 'No region found for lon=',str(lon),', lat=',str(lat)
+                    break #return
+                # Loop downward
+                for k in range(N):
+                    curr_temp = roms_temp[k,j,i]
+                    curr_salt = roms_salt[k,j,i]
+                    curr_volume = dV[k,j,i]
+                    # Get surface freezing point at this salinity
+                    curr_tfrz = curr_salt/(-18.48 + 18.48/1e3*curr_salt)
+                    # Figure out what water mass this is
+                    if curr_temp < curr_tfrz:
+                        # ISW
+                        wm_key = 0
+                    elif curr_salt < 34:
+                        # AASW
+                        wm_key = 4
+                    elif curr_temp > -1.5:
+                        # MCDW
+                        wm_key = 1
+                    elif curr_salt < 34.5:
+                        # LSSW
+                        wm_key = 3
+                    else:
+                        # HSSW
+                        wm_key = 2
+                    # Integrate volume for the right water mass and sector
+                    roms_vol_watermass[wm_key, sector] += curr_volume
+                    # Also integrate total Antarctica
+                    roms_vol_watermass[wm_key, -1] += curr_volume
+    # Find total volume of each sector by adding up the volume of each
+    # water mass
+    roms_vol_sectors = sum(roms_vol_watermass, axis=0)
+    # Calculate percentage of each water mass in each sector
+    roms_percent_watermass = zeros([num_watermasses, num_sectors])
+    for wm_key in range(num_watermasses):
+        for sector in range(num_sectors):
+            roms_percent_watermass[wm_key, sector] = roms_vol_watermass[wm_key, sector]/roms_vol_sectors[sector]*100                
+
+    print 'Processing low-res FESOM'
+    # Build mesh
+    elements_lr = fesom_grid(fesom_mesh_lr, circumpolar, cross_180)
+    id = Dataset(fesom_file_lr, 'r')
+    temp_nodes_lr = id.variables['temp'][0,:]
+    salt_nodes_lr = id.variables['salt'][0,:]
+    id.close()
+    fesom_vol_watermass_lr = zeros([num_watermasses, num_sectors])
+    for i in range(len(elements_lr)):
+        elm = elements_lr[i]
+        if elm.cavity:
+            lon = mean(elm.lon)
+            lat = mean(elm.lat)
+            if lon >= -85 and lon < -30 and lat < -74:
+                sector = 0
+            elif lon >= -30 and lon < 65:
+                sector = 1
+            elif lon >= 65 and lon < 76:
+                sector = 2
+            elif lon >= 76 and lon < 165 and lat >= -74:
+                sector = 3
+            elif (lon >= 155 and lon < 165 and lat < -74) or (lon >= 165) or (lon < -140):
+                sector = 4
+            elif (lon >= -140 and lon < -105) or (lon >= -105 and lon < -98 and lat < -73.1):
+                sector = 5
+            elif (lon >= -104 and lon < -98 and lat >= -73.1) or (lon >= -98 and lon < -66 and lat >= -75):
+                sector = 6
+            elif lon >= -66 and lon < -59 and lat >= -74:
+                sector = 7
+            else:
+                print 'No region found for lon=',str(lon),', lat=',str(lat)
+                break #return
+            # Get area of 2D element
+            area = elm.area()
+            nodes = [elm.nodes[0], elm.nodes[1], elm.nodes[2]]
+            # Loop downward
+            while True:
+                if nodes[0].below is None or nodes[1].below is None or nodes[2].below is None:
+                    # Reached the bottom
+                    break
+                # Calculate average temperature, salinity, and
+                # layer thickness for this 3D triangular prism
+                temp_vals = []
+                salt_vals = []
+                dz_vals = []
+                for n in range(3):
+                    temp_vals.append(temp_nodes_lr[nodes[n].id])
+                    salt_vals.append(salt_nodes_lr[nodes[n].id])
+                    temp_vals.append(temp_nodes_lr[nodes[n].below.id])
+                    salt_vals.append(salt_nodes_lr[nodes[n].below.id])
+                    dz_vals.append(abs(nodes[n].depth - nodes[n].below.depth))
+                    # Get ready for next iteration of loop
+                    nodes[n] = nodes[n].below
+                curr_temp = mean(array(temp_vals))
+                curr_salt = mean(array(salt_vals))
+                curr_volume = area*mean(array(dz_vals))
+                curr_tfrz = -0.0575*curr_salt + 1.7105e-3*sqrt(curr_salt**3) - 2.155e-4*curr_salt**2
+                if curr_temp < curr_tfrz:
+                    wm_key = 0
+                elif curr_salt < 34:
+                    wm_key = 4
+                elif curr_temp > -1.5:
+                    wm_key = 1
+                elif curr_salt < 34.5:
+                    wm_key = 3
+                else:
+                    wm_key = 2
+                fesom_vol_watermass_lr[wm_key, sector] += curr_volume
+                fesom_vol_watermass_lr[wm_key, -1] += curr_volume
+    fesom_vol_sectors_lr = sum(fesom_vol_watermass_lr, axis=0)
+    fesom_percent_watermass_lr = zeros([num_watermasses, num_sectors])
+    for wm_key in range(num_watermasses):
+        for sector in range(num_sectors):
+            fesom_percent_watermass_lr[wm_key, sector] = fesom_vol_watermass_lr[wm_key, sector]/fesom_vol_sectors_lr[sector]*100
+
+    print 'Processing high-res FESOM'
+    elements_hr = fesom_grid(fesom_mesh_hr, circumpolar, cross_180)
+    fesom_vol_watermass_hr = zeros([num_watermasses, num_sectors])
+    id = Dataset(fesom_file_hr, 'r')
+    temp_nodes_hr = id.variables['temp'][0,:]
+    salt_nodes_hr = id.variables['salt'][0,:]
+    id.close()
+    for i in range(len(elements_hr)):
+        elm = elements_hr[i]
+        if elm.cavity:
+            lon = mean(elm.lon)
+            lat = mean(elm.lat)
+            if lon >= -85 and lon < -30 and lat < -74:
+                sector = 0
+            elif lon >= -30 and lon < 65:
+                sector = 1
+            elif lon >= 65 and lon < 76:
+                sector = 2
+            elif lon >= 76 and lon < 165 and lat >= -74:
+                sector = 3
+            elif (lon >= 155 and lon < 165 and lat < -74) or (lon >= 165) or (lon < -140):
+                sector = 4
+            elif (lon >= -140 and lon < -105) or (lon >= -105 and lon < -98 and lat < -73.1):
+                sector = 5
+            elif (lon >= -104 and lon < -98 and lat >= -73.1) or (lon >= -98 and lon < -66 and lat >= -75):
+                sector = 6
+            elif lon >= -66 and lon < -59 and lat >= -74:
+                sector = 7
+            else:
+                print 'No region found for lon=',str(lon),', lat=',str(lat)
+                break #return
+            area = elm.area()
+            nodes = [elm.nodes[0], elm.nodes[1], elm.nodes[2]]
+            while True:
+                if nodes[0].below is None or nodes[1].below is None or nodes[2].below is None:
+                    break
+                temp_vals = []
+                salt_vals = []
+                dz_vals = []
+                for n in range(3):
+                    temp_vals.append(temp_nodes_hr[nodes[n].id])
+                    salt_vals.append(salt_nodes_hr[nodes[n].id])
+                    temp_vals.append(temp_nodes_hr[nodes[n].below.id])
+                    salt_vals.append(salt_nodes_hr[nodes[n].below.id])
+                    dz_vals.append(abs(nodes[n].depth - nodes[n].below.depth))
+                    nodes[n] = nodes[n].below
+                curr_temp = mean(array(temp_vals))
+                curr_salt = mean(array(salt_vals))
+                curr_volume = area*mean(array(dz_vals))
+                curr_tfrz = -0.0575*curr_salt + 1.7105e-3*sqrt(curr_salt**3) - 2.155e-4*curr_salt**2
+                if curr_temp < curr_tfrz:
+                    wm_key = 0
+                elif curr_salt < 34:
+                    wm_key = 4
+                elif curr_temp > -1.5:
+                    wm_key = 1
+                elif curr_salt < 34.5:
+                    wm_key = 3
+                else:
+                    wm_key = 2
+                fesom_vol_watermass_hr[wm_key, sector] += curr_volume
+                fesom_vol_watermass_hr[wm_key, -1] += curr_volume
+    fesom_vol_sectors_hr = sum(fesom_vol_watermass_hr, axis=0)
+    fesom_percent_watermass_hr = zeros([num_watermasses, num_sectors])
+    for wm_key in range(num_watermasses):
+        for sector in range(num_sectors):
+            fesom_percent_watermass_hr[wm_key, sector] = fesom_vol_watermass_hr[wm_key, sector]/fesom_vol_sectors_hr[sector]*100
+
+    print 'Plotting'
+    fig = figure(figsize=(12,6))
+    gs = GridSpec(3,3)
+    gs.update(left=0.15, right=0.98, bottom=0.2, top=0.88, wspace=0.1, hspace=0.28)
+    handles = []
+    for sector in range(num_sectors):
+        ax = subplot(gs[sector/3, sector%3])
+        lefts = 0
+        for wm_key in range(num_watermasses):
+            ax.barh(0, roms_percent_watermass[wm_key, sector], color=wm_colours[wm_key], left=lefts, align='center')
+            lefts += roms_percent_watermass[wm_key, sector]
+        lefts = 0
+        for wm_key in range(num_watermasses):
+            ax.barh(1, fesom_percent_watermass_lr[wm_key, sector], color=wm_colours[wm_key], left=lefts, align='center')
+            lefts += fesom_percent_watermass_lr[wm_key, sector]
+        lefts = 0
+        for wm_key in range(num_watermasses):
+            tmp = ax.barh(2, fesom_percent_watermass_hr[wm_key, sector], color=wm_colours[wm_key], left=lefts, align='center')
+            if sector == num_sectors-1:
+                handles.append(tmp)
+            lefts += fesom_percent_watermass_hr[wm_key, sector]
+        xlim([0, 100])
+        ax.invert_yaxis()
+        ax.set_yticks(range(3))
+        if sector % 3 == 0:
+            ax.set_yticklabels(('MetROMS', 'FESOM (low-res)', 'FESOM (high-res)'))
+        else:
+            ax.set_yticklabels(('','',''))
+        if sector >= num_sectors-3:
+            ax.set_xlabel('% volume')
+        else:
+            ax.set_xticklabels([])
+        ax.set_title(sector_names[sector])
+    legend(handles, wm_names, ncol=num_watermasses, bbox_to_anchor=(0.35,-0.4))
+    subplots_adjust(wspace=0.05, hspace=0.2)
+    suptitle('Water masses in ice shelf cavities (2002-2016 average)', fontsize=20)
+    fig.show()
+    fig.savefig('wm_barchart.png')
+
+
+# Command-line interface
+if __name__ == "__main__":
+
+    roms_grid = raw_input("Path to ROMS grid file: ")
+    roms_file = raw_input("Path to ROMS time-averaged temperature and salinity file, 2002-2016: ")
+    fesom_mesh_lr = raw_input("Path to FESOM low-res mesh directory: ")
+    fesom_file_lr = raw_input("Path to FESOM low-res time-averaged temperature and salinity file, 2002-2016: ")
+    fesom_mesh_hr = raw_input("Path to FESOM high-res mesh directory: ")
+    fesom_file_hr = raw_input("Path to FESOM high-res time-averaged temperature and salinity file, 2002-2016: ")
+    mip_watermass_barchart(roms_grid, roms_file, fesom_mesh_lr, fesom_mesh_hr, fesom_file_lr, fesom_file_hr)
+
+    
+    
+
+    
+    
+    
+    
diff --git a/timeseries_massloss_bellingshausen.py b/timeseries_massloss_bellingshausen.py
new file mode 100644
index 0000000..a9dac5c
--- /dev/null
+++ b/timeseries_massloss_bellingshausen.py
@@ -0,0 +1,118 @@
+from netCDF4 import Dataset
+from numpy import *
+from timeseries_massloss import calc_grid
+
+# Calculate timeseries of basal mass loss for the Wilkins, Stange, and
+# George VI ice shelves separately. Output to a log file.
+# Input:
+# directory = path to directory containing containing ocean averages files,
+#             assuming 5-day averages
+# start_index, end_index = integers containing range of files to process. For
+#                          example, start_index=1 and end_index=102 will
+#                          process files ocean_avg_0001.nc through
+#                          ocean_avg_0102.nc.
+# log_path = path to desired output log file
+def timeseries_massloss_bellingshausen (directory, start_index, end_index, log_path):
+
+    # Name of each ice shelf
+    names = ['Wilkins Ice Shelf', 'Stange Ice Shelf', 'George VI Ice Shelf']
+    # Limits on longitude and latitude for each ice shelf
+    # These depend on the source geometry, in this case RTopo 1.05
+    # Note there is one extra index at the end of each array; this is because
+    # the George VI region is awkwardly shaped and split into 2 (and hence why
+    # all 3 ice shelves are lumped together in timeseries_massloss.py)
+    lon_min = [-75, -79, -74.5, -69.5]
+    lon_max = [-69, -74.5, -67, -66]
+    lat_min = [-71.5, -73.8, -73.8, -72.6]
+    lat_max = [-69, -72.6, -72.6, -70]
+    # Density of ice in kg/m^3
+    rho_ice = 916
+
+    print 'Analysing grid'
+    # Using the first file, calculate dA (masked with ice shelf mask) and
+    # lon and lat coordinates
+    dA, lon, lat = calc_grid(directory + index_to_file(start_index))
+
+    print 'Setting up arrays'
+    # Set up array of masked area values for each ice shelf
+    dA_masked = ma.empty([len(names), size(dA,0), size(dA,1)])
+    for shelf in range(len(names)):
+        # Mask dA for the current ice shelf (note dA is already masked
+        # with the global ice shelf mask, so just restrict the indices
+        # to isolate the given ice shelf)
+        if shelf == len(names)-1:
+            # George VI region is split into two
+            dA_tmp = ma.masked_where(((lon < lon_min[shelf]) + (lon > lon_max[shelf]) + (lat < lat_min[shelf]) + (lat > lat_max[shelf]))*((lon < lon_min[shelf+1]) + (lon > lon_max[shelf+1]) + (lat < lat_min[shelf+1]) + (lat > lat_max[shelf+1])), dA)
+        else:
+            dA_tmp = ma.masked_where((lon < lon_min[shelf]) + (lon > lon_max[shelf]) + (lat < lat_min[shelf]) + (lat > lat_max[shelf]), dA)
+        dA_masked[shelf,:,:] = dA_tmp[:,:]
+    # Figure out how many time indices in total there are
+    total_time = 0
+    for index in range(start_index, end_index+1):
+        id = Dataset(directory + index_to_file(index), 'r')
+        total_time += id.variables['ocean_time'].size
+        id.close()
+    # Set up array of mass loss values for each ice shelf
+    massloss = empty([len(names), total_time])
+    # Also time values
+    time = empty(total_time)
+
+    t_posn = 0
+    for index in range(start_index, end_index+1):
+        file = directory + index_to_file(index)
+        print 'Processing ' + file
+        id = Dataset(file, 'r')
+        # Read melt rate and convert from m/s to m/y
+        ismr = id.variables['m'][:,:-15,1:-1]*365.25*24*60*60
+        # Read time and convert from seconds to years
+        curr_time = id.variables['ocean_time'][:]/(365.25*24*60*60)
+        id.close()
+        # Save the time values
+        num_time = size(curr_time)
+        time[t_posn:t_posn+num_time] = curr_time
+        for t in range(num_time):
+            for shelf in range(len(names)):
+                # Integrate ice shelf melt rate over area to get volume loss
+                volumeloss = sum(ismr[t,:,:]*dA_masked[shelf,:,:])
+                # Convert to mass loss in Gt/y
+                massloss[shelf, t_posn+t] = 1e-12*rho_ice*volumeloss
+        t_posn += num_time
+
+    print 'Saving results to log file'
+    f = open(log_path, 'w')
+    f.write('Time (years):\n')
+    for t in range(size(time)):
+        f.write(str(time[t]) + '\n')
+    for shelf in range(len(names)):
+        f.write(names[shelf] + ' Basal Mass Loss\n')
+        for t in range(size(time)):
+            f.write(str(massloss[shelf, t]) + '\n')
+    f.close()
+    
+
+# Given an integer, return the filename for the corresponding ocean averages
+# file. For example, index_to_file(1) = 'ocean_avg_0001.nc', and
+# index_to_file(95) = 'ocean_avg_0095.nc'.
+def index_to_file (index):
+
+    if index < 10:
+        return 'ocean_avg_000' + str(index) + '.nc'
+    elif index < 100:
+        return 'ocean_avg_00' + str(index) + '.nc'
+    elif index < 1000:
+        return 'ocean_avg_0' + str(index) + '.nc'
+    else:
+        return 'ocean_avg_' + str(index) + '.nc'
+
+
+# Command-line interface
+if __name__ == "__main__":
+
+    directory = raw_input("Path to ROMS output directory: ")
+    start_index = int(raw_input("Index of first ocean averages file (e.g. ocean_avg_0001.nc has index 1): "))
+    end_index = int(raw_input("Index of last ocean averages file: "))
+    log_path = raw_input("Path to desired logfile: ")
+    timeseries_massloss_bellingshausen(directory, start_index, end_index, log_path)
+
+
+