Removed northern boundary sponge layer and periodic boundary overlap …

…from plots and calculations

add_iceberg_melt.py

@@ -58,7 +58,7 @@ def add_iceberg_melt (file):
         melt_roms = melt_roms/rho_w*seconds_per_12h
         # Add to precipitation field for this month
         id = Dataset(file, 'a')
-        id.variables['rain'][12*year+month,:,:] = id.variables['rain'][12*year+month,:,:] - melt_roms
+        id.variables['rain'][month,:,:] = id.variables['rain'][month,:,:] - melt_roms
         id.close()
 
 

aice_animation.py

@@ -13,7 +13,7 @@
 # Directory containing CICE output files
 directory = '/short/y99/kaa561/roms_spinup_newest/cice/'
 # Number of time indices in each file
-num_ts = [144]
+num_ts = [504]
 # File number to start with for the animation (1-based)
 start_file = 1
 # Degrees to radians conversion factor
@@ -23,8 +23,8 @@
 
 # Read grid from the first file
 id = Dataset(directory + 'iceh' + str(start_file) + '.nc', 'r')
-lon = id.variables['TLON'][:,:]
-lat = id.variables['TLAT'][:,:]
+lon = id.variables['TLON'][:-15,:]
+lat = id.variables['TLAT'][:-15,:]
 id.close()
 
 # Calculate x and y coordinates for polar projection
@@ -53,7 +53,7 @@ def animate(i):
     id = Dataset(directory + 'iceh' + str(file_num) + '.nc', 'r')
     time_id = id.variables['time']
     time = num2date(time_id[i], units=time_id.units, calendar=time_id.calendar.lower()) 
-    data = id.variables['aice'][i,:,:]
+    data = id.variables['aice'][i,:-15,:]
     id.close()
     # Clear plot to save memory
     ax.collections = []
@@ -65,6 +65,6 @@ def animate(i):
     return img
 
 # Animate once every time index from start_file to the last file
-anim = FuncAnimation(fig, func=animate, frames=range(71,144)) #sum(array(num_ts[start_file-1:])))
+anim = FuncAnimation(fig, func=animate, frames=range(431,504)) #sum(array(num_ts[start_file-1:])))
 # Save as an mp4 with one frame per second
 anim.save('aice.mp4', fps=1)

avg_zeta.py

@@ -12,9 +12,9 @@ def avg_zeta (file_path):
     time = file.variables['ocean_time'][:]
     # Convert time from seconds to years
     time = time/(365*24*60*60)
-    lon = file.variables['lon_rho'][:,:]
-    lat = file.variables['lat_rho'][:,:]
-    mask = file.variables['mask_rho'][:,:]
+    lon = file.variables['lon_rho'][:-15,:-3]
+    lat = file.variables['lat_rho'][:-15,:-3]
+    mask = file.variables['mask_rho'][:-15,:-3]
     avg_zeta = []
 
     # Calculate dx and dy in another script
@@ -26,7 +26,7 @@ def avg_zeta (file_path):
     for l in range(size(time)):
         print 'Processing timestep ' + str(l+1) + ' of ' + str(size(time))
         # Read zeta at this timestep
-        zeta = file.variables['zeta'][l,:,:]
+        zeta = file.variables['zeta'][l,:-15,:-3]
         # Calculate area-weighted average
         avg_zeta.append(sum(zeta*dA)/sum(dA))
 

bwsalt_plot.py

@@ -0,0 +1,85 @@
+from netCDF4 import Dataset
+from numpy import *
+from matplotlib.pyplot import *
+
+# Create a circumpolar plot of bottom water salinity, averaged over the last
+# year of simulation.
+# Input:
+# file_path = path to ocean averages file containing at least one year of
+#             5-day averages
+# save = optional boolean to save the figure to a file, rather than display it
+#        on screen
+# fig_name = if save=True, filename for figure
+def bwsalt_plot (file_path, save=False, fig_name=None):
+
+    # Degrees to radians conversion factor
+    deg2rad = pi/180
+    # Centre of missing circle in grid
+    lon_c = 50
+    lat_c = -83
+    # Radius of missing circle (play around with this until it works)
+    radius = 10.1
+    # Bounds on colour scale
+    min_scale = 34
+    max_scale = 35
+
+    # Read the grid
+    id = Dataset(file_path, 'r')
+    lon = id.variables['lon_rho'][:-15,:-2]
+    lat = id.variables['lat_rho'][:-15,:-2]
+    # Read the last year of bottom water salinity (assume 5-day averages here)
+    # and average over time
+    bwsalt = mean(id.variables['salt'][-73:,0,:-15,:-2], axis=0)
+    id.close()
+
+    # Convert grid to spherical coordinates
+    x = -(lat+90)*cos(lon*deg2rad+pi/2)
+    y = (lat+90)*sin(lon*deg2rad+pi/2)
+    # Find centre in spherical coordinates
+    x_c = -(lat_c+90)*cos(lon_c*deg2rad+pi/2)
+    y_c = (lat_c+90)*sin(lon_c*deg2rad+pi/2)
+    # Build a regular x-y grid and select the missing circle
+    x_reg, y_reg = meshgrid(linspace(amin(x), amax(x), num=1000), linspace(amin(y), amax(y), num=1000))
+    land_circle = zeros(shape(x_reg))
+    land_circle = ma.masked_where(sqrt((x_reg-x_c)**2 + (y_reg-y_c)**2) > radius, land_circle)
+
+    # Set colour bounds
+    lev = linspace(min_scale, max_scale, num=50)
+
+    # Plot
+    fig = figure(figsize=(16,12))
+    ax = fig.add_subplot(1,1,1,aspect='equal')
+    fig.patch.set_facecolor('white')
+    # First shade everything in grey
+    contourf(x, y, ones(shape(bwsalt)), 1, colors=(('0.6', '0.6', '0.6')))
+    # Fill in the missing circle
+    contourf(x_reg, y_reg, land_circle, 1, colors=(('0.6', '0.6', '0.6')))
+    # Now shade the salinity (land mask will remain grey)
+    contourf(x, y, bwsalt, lev, cmap='jet', extend='both')
+    cbar = colorbar(ticks=arange(min_scale, max_scale+0.2, 0.2))
+    cbar.ax.tick_params(labelsize=20)
+    title(r'Bottom water salintiy (psu), annual average', fontsize=30)
+    axis('off')
+
+    # Finished
+    if save:
+        fig.savefig(fig_name)
+    else:
+        fig.show()
+
+
+# Command-line interface
+if __name__ == '__main__':
+
+    file_path = raw_input("Path to ocean averages file, containing at least 1 year of 5-day averages: ")
+    action = raw_input("Save figure (s) or display in window (d)? ")
+    if action == 's':
+        save = True
+        fig_name = raw_input("File name for figure: ")
+    else:
+        save = False
+        fig_name = None
+    # Make the plot
+    bwsalt_plot(file_path, save, fig_name)
+
+

bwtemp_plot.py

@@ -25,11 +25,11 @@ def bwtemp_plot (file_path, save=False, fig_name=None):
 
     # Read the grid
     id = Dataset(file_path, 'r')
-    lon = id.variables['lon_rho'][:,:]
-    lat = id.variables['lat_rho'][:,:]
+    lon = id.variables['lon_rho'][:-15,:-2]
+    lat = id.variables['lat_rho'][:-15,:-2]
     # Read the last year of bottom water temp (assume 5-day averages here) and
     # average over time
-    bwtemp = mean(id.variables['temp'][-73:,0,:,:], axis=0)
+    bwtemp = mean(id.variables['temp'][-73:,0,:-15,:-2], axis=0)
     id.close()
 
     # Convert grid to spherical coordinates

circumpolar_cice_plot.py

@@ -19,7 +19,7 @@ def circumpolar_cice_plot (file_path, var_name, tstep, colour_bounds=None, save=
 
     # Read the variable and figure out which grid it's on
     id = Dataset(file_path, 'r')
-    data = id.variables[var_name][tstep-1,:,:]
+    data = id.variables[var_name][tstep-1,:-15,:]
     if var_name == 'aice':
         units = 'fraction'
     else:
@@ -39,8 +39,8 @@ def circumpolar_cice_plot (file_path, var_name, tstep, colour_bounds=None, save=
         return
 
     # Read the correct lat and lon for this grid
-    lon = id.variables[lon_name][:,:]
-    lat = id.variables[lat_name][:,:]
+    lon = id.variables[lon_name][:-15,:]
+    lat = id.variables[lat_name][:-15,:]
     id.close()
 
     # Convert to spherical coordinates

circumpolar_plot.py

@@ -36,11 +36,11 @@ def circumpolar_plot (file_path, var_name, tstep, depth_key, depth, depth_bounds
     id = Dataset(file_path, 'r')
     if len(id.variables[var_name].shape) == 4:
         # 3D variable; will have to choose depth later
-        data_full = id.variables[var_name][tstep-1,:,:,:]
+        data_full = id.variables[var_name][tstep-1,:,:-15,:]
         choose_depth = True
     elif len(id.variables[var_name].shape) == 3:
         # 2D variable
-        data = id.variables[var_name][tstep-1,:,:]
+        data = id.variables[var_name][tstep-1,:-15,:]
         choose_depth = False
     if var_name == 'salt':
         units = 'psu'
@@ -72,8 +72,8 @@ def circumpolar_plot (file_path, var_name, tstep, depth_key, depth, depth_bounds
         return
 
     # Read grid variables
-    h = id.variables['h'][:,:]
-    zice = id.variables['zice'][:,:]
+    h = id.variables['h'][:-15,:]
+    zice = id.variables['zice'][:-15,:]
     # h and zice are on the rho-grid; interpolate if necessary
     if grid_name == 'u':
         h = 0.5*(h[:,0:-1] + h[:,1:])
@@ -82,10 +82,30 @@ def circumpolar_plot (file_path, var_name, tstep, depth_key, depth, depth_bounds
         h = 0.5*(h[0:-1,:] + h[1:,:])
         zice = 0.5*(zice[0:-1,:] + zice[1:,:])
     # Read the correct lat and lon for this grid
-    lon = id.variables[lon_name][:,:]
-    lat = id.variables[lat_name][:,:]
+    lon = id.variables[lon_name][:-15,:]
+    lat = id.variables[lat_name][:-15,:]
     id.close()
 
+    # Throw away the overlapping periodic boundary
+    if grid_name == 'u':
+        if choose_depth:
+            data_full = data_full[:,:,:-1]
+        else:
+            data = data[:,:-1]
+        lon = lon[:,:-1]
+        lat = lat[:,:-1]
+        h = h[:,:-1]
+        zice = zice[:,:-1]
+    else:
+        if choose_depth:
+            data_full = data_full[:,:,:-2]
+        else:
+            data = data[:,:-2]
+        lon = lon[:,:-2]
+        lat = lat[:,:-2]
+        h = h[:,:-2]
+        zice = zice[:,:-2]
+
     # Convert to spherical coordinates
     x = -(lat+90)*cos(lon*deg2rad+pi/2)
     y = (lat+90)*sin(lon*deg2rad+pi/2)

dpt_2d.py

@@ -31,11 +31,11 @@ def dpt_2d (file_path):
     print 'Reading grid'
     # Read grid variables
     id = Dataset(file_path, 'r')
-    h = id.variables['h'][:,:]
-    zice = id.variables['zice'][:,:]
-    lon = id.variables['lon_rho'][:,:]
-    lat = id.variables['lat_rho'][:,:]
-    mask = id.variables['mask_rho'][:,:]
+    h = id.variables['h'][:-15,:-3]
+    zice = id.variables['zice'][:-15,:-3]
+    lon = id.variables['lon_rho'][:-15,:-3]
+    lat = id.variables['lat_rho'][:-15,:-3]
+    mask = id.variables['mask_rho'][:-15,:-3]
 
     # Interpolate latitude to the edges of each cell
     s_bdry = lat[0,:]
@@ -67,11 +67,13 @@ def dpt_2d (file_path):
 
         print 'Processing timestep ' + str(t+1) + ' of '+str(size(time))
         # Read ubar and interpolate onto the rho-grid
-        ubar = id.variables['ubar'][t,:,:]
+        ubar = id.variables['ubar'][t,:-15,:]
         w_bdry_ubar = 0.5*(ubar[:,0] + ubar[:,-1])
         middle_ubar = 0.5*(ubar[:,0:-1] + ubar[:,1:])
         e_bdry_ubar = w_bdry_ubar[:]
         ubar_rho = ma.concatenate((w_bdry_ubar[:,None], middle_ubar, e_bdry_ubar[:,None]), axis=1)
+        # Throw away the overlapping periodic boundary
+        ubar_rho = ubar_rho[:,:-3]
         # Trim to Drake Passage bounds
         ubar_rho_DP = ubar_rho[j_min:j_max,i_DP]
         # Calculate transport and convert to Sv

dpt_2d_int.py

@@ -30,11 +30,11 @@ def dpt_2d_int (file_path):
     print 'Reading grid'
     # Read grid variables
     id = Dataset(file_path, 'r')
-    h = id.variables['h'][:,:]
-    zice = id.variables['zice'][:,:]
-    lon = id.variables['lon_rho'][:,:]
-    lat = id.variables['lat_rho'][:,:]
-    mask = id.variables['mask_rho'][:,:]
+    h = id.variables['h'][:-15,:-3]
+    zice = id.variables['zice'][:-15,:-3]
+    lon = id.variables['lon_rho'][:-15,:-3]
+    lat = id.variables['lat_rho'][:-15,:-3]
+    mask = id.variables['mask_rho'][:-15,:-3]
 
     # Interpolate latitude to the edges of each cell
     s_bdry = lat[0,:]
@@ -66,11 +66,13 @@ def dpt_2d_int (file_path):
 
         print 'Processing timestep ' + str(t+1) + ' of '+str(size(time))
         # Read ubar and interpolate onto the rho-grid
-        ubar = id.variables['ubar'][t,:,:]
+        ubar = id.variables['ubar'][t,:-15,:]
         w_bdry_ubar = 0.5*(ubar[:,0] + ubar[:,-1])
         middle_ubar = 0.5*(ubar[:,0:-1] + ubar[:,1:])
         e_bdry_ubar = w_bdry_ubar[:]
         ubar_rho = ma.concatenate((w_bdry_ubar[:,None], middle_ubar, e_bdry_ubar[:,None]), axis=1)
+        # Throw away the overlapping periodic boundary
+        ubar_rho = ubar_rho[:,:-3]
         # Trim to Drake Passage bounds
         ubar_rho_DP = ubar_rho[j_min:j_max,i_DP]
         # Calculate transport and convert to Sv

dpt_timeseries.py

@@ -37,11 +37,11 @@ def dpt_timeseries (file_path):
     print 'Reading grid'
     # Read grid variables
     id = Dataset(file_path, 'r')
-    h = id.variables['h'][:,:]
-    zice = id.variables['zice'][:,:]
-    lon = id.variables['lon_rho'][:,:]
-    lat = id.variables['lat_rho'][:,:]
-    mask = id.variables['mask_rho'][:,:]
+    h = id.variables['h'][:-15,:-3]
+    zice = id.variables['zice'][:-15,:-3]
+    lon = id.variables['lon_rho'][:-15,:-3]
+    lat = id.variables['lat_rho'][:-15,:-3]
+    mask = id.variables['mask_rho'][:-15,:-3]
 
 
     # Interpolate latitude to the edges of each cell
@@ -73,11 +73,13 @@ def dpt_timeseries (file_path):
 
         print 'Processing timestep ' + str(t+1) + ' of '+str(size(time))
         # Read ubar and interpolate onto the rho-grid
-        ubar = id.variables['ubar'][t,:,:]
+        ubar = id.variables['ubar'][t,:-15,:]
         w_bdry_ubar = 0.5*(ubar[:,0] + ubar[:,-1])
         middle_ubar = 0.5*(ubar[:,0:-1] + ubar[:,1:])
         e_bdry_ubar = w_bdry_ubar[:]
         ubar_rho = ma.concatenate((w_bdry_ubar[:,None], middle_ubar, e_bdry_ubar[:,None]), axis=1)
+        # Throw away the overlapping periodic boundary
+        ubar_rho = ubar_rho[:,:-3]
         # Trim to Drake Passage bounds
         ubar_rho_DP = ubar_rho[j_min:j_max,i_DP]
         # Integrate transport and convert to Sv

file_guide.txt

@@ -585,6 +585,15 @@ bwtemp_plot.py: Creates a circumpolar plot of bottom water temperature, averaged
 			figure (and if so, what filename) or display it on
 			screen.
 
+bwsalt_plot.py: Creates a circumpolar plot of bottom water salinity, averaged
+                over the last year of simulation.
+		To run: Open python or ipython and type "run bwsalt_plot.py".
+		        The script will prompt you for the path to a ROMS
+			output file containing at least one year of 5-day
+			averages, and ask you whether you want to save the
+			figure (and if so, what filename) or display it on
+			screen.
+
 ismr_seasonal_cycle.py: Make a map of the magnitude of the seasonal cycle (over
                         the last year of simulation) in area-averaged melt rates
                         for each ice shelf that is over 5,000 km^2 in Rignot

h_circumpolar.py

@@ -19,10 +19,10 @@ def h_circumpolar (grid_path, fig_name):
 
     # Read data
     id = Dataset(grid_path, 'r')
-    data = id.variables['h'][:,:]
-    lon = id.variables['lon_rho'][:,:]
-    lat = id.variables['lat_rho'][:,:]
-    mask = id.variables['mask_rho'][:,:]
+    data = id.variables['h'][:-15,:-2]
+    lon = id.variables['lon_rho'][:-15,:-2]
+    lat = id.variables['lat_rho'][:-15,:-2]
+    mask = id.variables['mask_rho'][:-15,:-2]
     id.close()
 
     # Mask with land mask

ice2ocn_fwflux.py

@@ -16,9 +16,9 @@ def ice2ocn_fwflux (file_path, tstep, fig_name):
 
     # Read data
     id = Dataset(file_path, 'r')
-    data = id.variables['fresh_ai'][tstep-1,:,:] - id.variables['fsalt_ai'][tstep-1,:,:]/1000.*60.*60.*24.*100.
-    lon = id.variables['TLON'][:,:]
-    lat = id.variables['TLAT'][:,:]
+    data = id.variables['fresh_ai'][tstep-1,:-15,:] - id.variables['fsalt_ai'][tstep-1,:-15,:]/1000.*60.*60.*24.*100.
+    lon = id.variables['TLON'][:-15,:]
+    lat = id.variables['TLAT'][:-15,:]
     id.close()
 
     # Convert to spherical coordinates
@@ -35,8 +35,8 @@ def ice2ocn_fwflux (file_path, tstep, fig_name):
     title('CICE-to-ROMS net freshwater flux (cm/day)', fontsize=30)
     axis('off')
 
-    #savefig(fig_name)
-    show()
+    savefig(fig_name)
+    #show()
 
 
 # Command-line interface