From be2997dbd4d14e67c68cd18783ab2a95b2593b84 Mon Sep 17 00:00:00 2001 From: Kaitlin Alexander Date: Fri, 5 Aug 2016 13:43:27 +1000 Subject: [PATCH] Removed northern boundary sponge layer and periodic boundary overlap from plots and calculations --- add_iceberg_melt.py | 2 +- aice_animation.py | 10 +-- avg_zeta.py | 8 +- bwsalt_plot.py | 85 +++++++++++++++++++ bwtemp_plot.py | 6 +- circumpolar_cice_plot.py | 6 +- circumpolar_plot.py | 32 ++++++-- dpt_2d.py | 14 ++-- dpt_2d_int.py | 14 ++-- dpt_timeseries.py | 14 ++-- file_guide.txt | 9 ++ h_circumpolar.py | 8 +- ice2ocn_fwflux.py | 10 +-- ini_sss_circumpolar.py | 9 +- ini_sst_circumpolar.py | 12 +-- ismr_map.py | 14 ++-- ismr_plot.py | 16 ++-- ismr_seasonal_cycle.py | 173 --------------------------------------- massloss.py | 8 +- massloss_map.py | 12 +-- max_vel.py | 4 +- nsidc_aice_monthly.py | 10 +-- nsidc_aice_seasonal.py | 22 ++--- other_models_data | 43 ++++++++++ overturning_plot.py | 15 ++-- sose_roms_seasonal.py | 14 ++-- spinup_plots.py | 40 +++++---- uv_vectorplot.py | 21 +++-- zice.py | 9 +- zice_circumpolar.py | 12 +-- zonal_plot.py | 28 +++++-- 31 files changed, 349 insertions(+), 331 deletions(-) create mode 100644 bwsalt_plot.py delete mode 100644 ismr_seasonal_cycle.py create mode 100644 other_models_data diff --git a/add_iceberg_melt.py b/add_iceberg_melt.py index c42ede0..57b9b7c 100644 --- a/add_iceberg_melt.py +++ b/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() diff --git a/aice_animation.py b/aice_animation.py index fc9876b..a3572be 100644 --- a/aice_animation.py +++ b/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) diff --git a/avg_zeta.py b/avg_zeta.py index ca104e4..3f56972 100644 --- a/avg_zeta.py +++ b/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)) diff --git a/bwsalt_plot.py b/bwsalt_plot.py new file mode 100644 index 0000000..dd14452 --- /dev/null +++ b/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) + + diff --git a/bwtemp_plot.py b/bwtemp_plot.py index 1718ab2..520ba9e 100644 --- a/bwtemp_plot.py +++ b/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 diff --git a/circumpolar_cice_plot.py b/circumpolar_cice_plot.py index 8a49576..b0e3e50 100644 --- a/circumpolar_cice_plot.py +++ b/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 diff --git a/circumpolar_plot.py b/circumpolar_plot.py index 4d86263..41ff2f9 100644 --- a/circumpolar_plot.py +++ b/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) diff --git a/dpt_2d.py b/dpt_2d.py index aa08173..f056cea 100644 --- a/dpt_2d.py +++ b/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 diff --git a/dpt_2d_int.py b/dpt_2d_int.py index 8c6ae88..655f3e8 100644 --- a/dpt_2d_int.py +++ b/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 diff --git a/dpt_timeseries.py b/dpt_timeseries.py index 9aef2a1..76aab92 100644 --- a/dpt_timeseries.py +++ b/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 diff --git a/file_guide.txt b/file_guide.txt index 0998937..f180068 100644 --- a/file_guide.txt +++ b/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 diff --git a/h_circumpolar.py b/h_circumpolar.py index dc20860..7015c34 100644 --- a/h_circumpolar.py +++ b/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 diff --git a/ice2ocn_fwflux.py b/ice2ocn_fwflux.py index 74a5d7d..1657d9c 100644 --- a/ice2ocn_fwflux.py +++ b/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 diff --git a/ini_sss_circumpolar.py b/ini_sss_circumpolar.py index cdc700a..b4b1fcf 100644 --- a/ini_sss_circumpolar.py +++ b/ini_sss_circumpolar.py @@ -22,14 +22,14 @@ def ini_sss_circumpolar (grid_path, file_path, fig_name): # Read surface temps id = Dataset(file_path, 'r') - data = id.variables['salt'][0,-1,:,:] + data = id.variables['salt'][0,-1,:-15,:-2] id.close() # Read lat and lon id = Dataset(grid_path, 'r') - lon = id.variables['lon_rho'][:,:] - lat = id.variables['lat_rho'][:,:] - mask = id.variables['mask_rho'][:,:] + 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 @@ -50,6 +50,7 @@ def ini_sss_circumpolar (grid_path, file_path, fig_name): axis('off') savefig(fig_name) + #show() # Command-line interface diff --git a/ini_sst_circumpolar.py b/ini_sst_circumpolar.py index b89ef94..c3b52fc 100644 --- a/ini_sst_circumpolar.py +++ b/ini_sst_circumpolar.py @@ -22,14 +22,14 @@ def ini_sst_circumpolar (grid_path, file_path, fig_name): # Read surface temps id = Dataset(file_path, 'r') - data = id.variables['temp'][0,-1,:,:] + data = id.variables['temp'][0,-1,:-15,:-2] id.close() # Read lat and lon id = Dataset(grid_path, 'r') - lon = id.variables['lon_rho'][:,:] - lat = id.variables['lat_rho'][:,:] - mask = id.variables['mask_rho'][:,:] + 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 @@ -49,8 +49,8 @@ def ini_sst_circumpolar (grid_path, file_path, fig_name): title(r'Initial SST ($^{\circ}$C)', fontsize=30) axis('off') - #savefig(fig_name) - show() + savefig(fig_name) + #show() # Command-line interface diff --git a/ismr_map.py b/ismr_map.py index dcb612f..c574be9 100644 --- a/ismr_map.py +++ b/ismr_map.py @@ -23,7 +23,7 @@ def ismr_map (grid_path, log_path, save=False, fig_name=None): lat_max = [-69.37, -66.13, -69.5, -74.67, -71.67, -74.17, -74.67, -73.67, -73, -75.17, -76.41, -66.67, -66.5, -64.83, -66.17, -68.33, -68.67, -68.33, -69.33, -69.83, -69.33, -71.5, -77.77, -77] # Area of each ice shelf in m^2 (printed to screen during massloss.py, # update if the grid changes) - area = [12754001970.4, 52008964915.9, 47287926558.6, 353435138251.0, 31290573250.5, 5162051654.52, 3990382861.08, 4680996769.75, 32446806852.2, 7694313052.38, 13537287121.0, 4918446447.87, 6482036686.01, 30521756982.6, 15158334399.6, 64359735004.9, 4575785549.65, 45327465354.5, 8110511960.62, 7088165282.99, 55762463981.5, 68006982027.4, 429252991746.0] + area = [12754001970.4, 52008964915.9, 47287926558.6, 353435138251.0, 31290573250.5, 5162051654.52, 3990382861.08, 4680996769.75, 32446806852.2, 7694313052.38, 13537287121.0, 4918446447.87, 6482036686.01, 30521756982.6, 15158334399.6, 64359735004.9, 4575785549.65, 45327465354.5, 8110511960.62, 7088165282.99, 54898163328.1, 68006982027.4, 429252991746.0] # Observed melt rate (Rignot 2013) and uncertainty for each ice shelf, in Gt/y obs_ismr = [0.1, 0.4, 3.1, 0.3, 1.7, 16.2, 17.7, 7.8, 4.3, 0.6, 1.5, 1.4, 7.7, 2.8, 1.7, 0.6, -0.4, 0.4, 0.7, 0.5, 0.5, 0.1, 0.1] obs_ismr_error = [0.6, 1, 0.8, 0.1, 0.6, 1, 1, 0.6, 0.4, 0.3, 0.3, 0.6, 0.7, 0.6, 0.7, 0.4, 0.6, 0.4, 0.2, 0.2, 0.2, 0.2, 0.1] @@ -80,11 +80,11 @@ def ismr_map (grid_path, log_path, save=False, fig_name=None): # Read the grid id = Dataset(grid_path, 'r') - lon = id.variables['lon_rho'][:,:] - lat = id.variables['lat_rho'][:,:] - mask_rho = id.variables['mask_rho'][:,:] - mask_zice = id.variables['mask_zice'][:,:] - zice = id.variables['zice'][:,:] + lon = id.variables['lon_rho'][:-15,:-2] + lat = id.variables['lat_rho'][:-15,:-2] + mask_rho = id.variables['mask_rho'][:-15,:-2] + mask_zice = id.variables['mask_zice'][:-15,:-2] + zice = id.variables['zice'][:-15,:-2] id.close() # Make sure longitude goes from -180 to 180, not 0 to 360 @@ -140,7 +140,7 @@ def ismr_map (grid_path, log_path, save=False, fig_name=None): # Determine bounds on colour scale max_val = amax(abs(error)) lev = linspace(-max_val, max_val, num=40) - lev = linspace(-200, 200, num=40) + lev = linspace(-100, 100, num=40) # Space ticks on colorbar 25% apart max_tick = floor(max_val/25)*25 diff --git a/ismr_plot.py b/ismr_plot.py index fdca12a..57086bc 100644 --- a/ismr_plot.py +++ b/ismr_plot.py @@ -29,25 +29,25 @@ def ismr_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'][:,:] - mask_rho = id.variables['mask_rho'][:,:] - zice = id.variables['zice'][:,:] + lon = id.variables['lon_rho'][:-15,:-2] + lat = id.variables['lat_rho'][:-15,:-2] + mask_rho = id.variables['mask_rho'][:-15,:-2] + zice = id.variables['zice'][:-15,:-2] # Read the last year of ice shelf melt rates (assume 5-day averages here), # average over time, and convert from m/s to m/y - ismr = mean(id.variables['m'][-73:,:,:], axis=0)*60*60*24*365.25 + ismr = mean(id.variables['m'][-73:,:-15,:-2], axis=0)*60*60*24*365.25 id.close() # Mask the open ocean and land out of the melt rates ismr = ma.masked_where(zice==0, ismr) # Set colour map # Values for change points - cmap_vals = array([-0.5, 0, 1, 2, 5, 8]) + cmap_vals = array([-0.1, 0, 1, 2, 5, 8]) # Colours for change points # (blue, white, yellow-orange, red-orange, dark red, purple) cmap_colors = [(0.26, 0.45, 0.86), (1, 1, 1), (1, 0.9, 0.4), (0.99, 0.59, 0.18), (0.5, 0.0, 0.08), (0.96, 0.17, 0.89)] # Map to 0-1 - cmap_vals_norm = (cmap_vals + 0.5)/(8 + 0.5) + cmap_vals_norm = (cmap_vals + 0.1)/(8 + 0.1) # Combine into a list cmap_list = [] for i in range(size(cmap_vals)): @@ -55,7 +55,7 @@ def ismr_plot (file_path, save=False, fig_name=None): # Make colour map mf_cmap = LinearSegmentedColormap.from_list('melt_freeze', cmap_list) # Set levels - lev = linspace(-0.5, 8, num=100) + lev = linspace(-0.1, 8, num=100) # Get land/zice mask open_ocn = copy(mask_rho) diff --git a/ismr_seasonal_cycle.py b/ismr_seasonal_cycle.py deleted file mode 100644 index 1f3ffbd..0000000 --- a/ismr_seasonal_cycle.py +++ /dev/null @@ -1,173 +0,0 @@ -from netCDF4 import Dataset -from numpy import * -from matplotlib.pyplot import * -from cartesian_grid_2d import * - -# 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 et al., 2013. -# Input: -# grid_path = path to ROMS grid file -# log_path = path to logfile created using massloss.py -# 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 ismr_seasonal_cycle (grid_path, log_path, save=False, fig_name=None): - - # 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 Ross region crosses the line 180W and therefore is split into two - lon_min = [-62.67, -65.5, -79.17, -85, -104.17, -102.5, -108.33, -114.5, -135.67, -149.17, -155, 144, 115, 94.17, 80.83, 65, 33.83, 19, 12.9, 9.33, -10.05, -28.33, -180, 158.33] - lon_max = [-59.33, -60, -66.67, -28.33, -88.83, -99.17, -103.33, -111.5, -114.33, -140, -145, 146.62, 123.33, 102.5, 89.17, 75, 37.67, 33.33, 16.17, 12.88, 7.6, -10.33, -146.67, 180] - lat_min = [-73.03, -69.35, -74.17, -83.5, -73.28, -75.5, -75.5, -75.33, -74.9, -76.42, -78, -67.83, -67.17, -66.67, -67.83, -73.67, -69.83, -71.67, -70.5, -70.75, -71.83, -76.33, -85, -84.5] - lat_max = [-69.37, -66.13, -69.5, -74.67, -71.67, -74.17, -74.67, -73.67, -73, -75.17, -76.41, -66.67, -66.5, -64.83, -66.17, -68.33, -68.67, -68.33, -69.33, -69.83, -69.33, -71.5, -77.77, -77] - num_iceshelves = len(lon_min)-1 - - # Density of ice in kg/m^3 - rho_ice = 916 - # Degrees to radians conversion factor - deg2rad = pi/180 - # Northern boundary 63S for plot - nbdry = -63+90 - # 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 - - # Read log file - time = [] - f = open(log_path, 'r') - # Skip the 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 - # Set up array for mass loss values at each ice shelf - massloss_ts = empty([num_iceshelves, len(time)]) - index = 0 - # Loop over ice shelves - while index < num_iceshelves: - t = 0 - for line in f: - try: - massloss_ts[index, t] = float(line) - t += 1 - except(ValueError): - # Reached the header for the next variable - break - index +=1 - - # Find the time indices we care about: last year of simulation - time = array(time) - min_t = nonzero(time >= time[-1]-1)[0][0] - max_t = size(time) - # Find the magnitude of the seasonal cycle in mass loss for each ice shelf - massloss_cycle = empty(num_iceshelves) - for index in range(num_iceshelves): - massloss_cycle[index] = amax(massloss_ts[index, min_t:max_t]) - amin(massloss_ts[index, min_t:max_t]) - - # Read the grid - id = Dataset(grid_path, 'r') - lon = id.variables['lon_rho'][:,:] - lat = id.variables['lat_rho'][:,:] - mask_rho = id.variables['mask_rho'][:,:] - mask_zice = id.variables['mask_zice'][:,:] - id.close() - - # Calculate dA and mask with zice - dx, dy = cartesian_grid_2d(lon, lat) - dA = ma.masked_where(mask_zice==0, dx*dy) - - # Make sure longitude goes from -180 to 180, not 0 to 360 - index = lon > 180 - lon[index] = lon[index] - 360 - - # Get land/zice mask - open_ocn = copy(mask_rho) - open_ocn[mask_zice==1] = 0 - land_zice = ma.masked_where(open_ocn==1, open_ocn) - - # Initialise a field of ice shelf melt rate seasonal cycle - ismr_cycle = zeros(shape(lon)) - # Mask out land and open ocean points - ismr_cycle = ma.masked_where(mask_zice==0, ismr_cycle) - - # Loop over ice shelves - for index in range(num_iceshelves): - # Find the grid cells within this ice shelf - if index == num_iceshelves-1: - # Ross region is split into two - region = (lon >= lon_min[index])*(lon <= lon_max[index])*(lat >= lat_min[index])*(lat <= lat_max[index])*(mask_zice == 1) + (lon >= lon_min[index+1])*(lon <= lon_max[index+1])*(lat >= lat_min[index+1])*(lat <= lat_max[index+1])*(mask_zice == 1) - else: - region = (lon >= lon_min[index])*(lon <= lon_max[index])*(lat >= lat_min[index])*(lat <= lat_max[index])*(mask_zice == 1) - # Calculate the conversion factor from mass loss to area-averaged - # melt rate for this ice shelf - factor = 1e12/(rho_ice*sum(dA[region])) - # Get seasonal cycle in melt rates and modify the field for this region - ismr_cycle[region] = massloss_cycle[index]*factor - - # 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(-nbdry, nbdry, num=1000), linspace(-nbdry, nbdry, 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 up colour levels - lev = linspace(0, amax(ismr_cycle), num=50) - print amax(ismr_cycle) - - # Set up plot - fig = figure(figsize=(16,12)) - ax = fig.add_subplot(1,1,1, aspect='equal') - fig.patch.set_facecolor('white') - # First shade land and zice in grey (include zice so there are no white - # patches near the grounding line where contours meet) - contourf(x, y, land_zice, 1, colors=(('0.6', '0.6', '0.6'))) - # Fill in the missing cicle - contourf(x_reg, y_reg, land_circle, 1, colors=(('0.6', '0.6', '0.6'))) - # Now shade the error in mass loss - contourf(x, y, ismr_cycle, lev, cmap='jet') - cbar = colorbar() - cbar.ax.tick_params(labelsize=20) - xlim([-nbdry, nbdry]) - ylim([-nbdry, nbdry]) - title('Seasonal cycle in area-averaged ice shelf melt rate (m/y)', fontsize=30) - axis('off') - - # Finished - if save: - fig.savefig(fig_name) - else: - fig.show() - - -# Command-line interface -if __name__ == "__main__": - - grid_path = raw_input("Path to grid file: ") - log_path = raw_input("Path to mass loss logfile: ") - action = raw_input("Save figure (s) or display in window (d)? ") - if action == 's': - save = True - fig_name = raw_input("File name for figure: ") - elif action == 'd': - save = False - fig_name = None - # Make the plot - ismr_seasonal_cycle(grid_path, log_path, save, fig_name) - - - - - - diff --git a/massloss.py b/massloss.py index 7bbeef4..ded151e 100644 --- a/massloss.py +++ b/massloss.py @@ -95,7 +95,7 @@ def massloss (file_path, log_path): # Read ice shelf melt rate, converting to float128 to prevent overflow # during integration - ismr = ma.asarray(id.variables['m'][t-start_t,:,:], dtype=float128) + ismr = ma.asarray(id.variables['m'][t-start_t,:-15,:-3], dtype=float128) # Convert from m/s to m/y ismr = ismr*365.25*24*60*60 @@ -199,9 +199,9 @@ def calc_grid (file_path): # Read grid variables id = Dataset(file_path, 'r') - lon = id.variables['lon_rho'][:,:] - lat = id.variables['lat_rho'][:,:] - zice = id.variables['zice'][:,:] + lon = id.variables['lon_rho'][:-15,:-3] + lat = id.variables['lat_rho'][:-15,:-3] + zice = id.variables['zice'][:-15,:-3] id.close() # Calculate dx and dy in another script diff --git a/massloss_map.py b/massloss_map.py index badc19c..fe849cb 100644 --- a/massloss_map.py +++ b/massloss_map.py @@ -74,11 +74,11 @@ def massloss_map (grid_path, log_path, save=False, fig_name=None): # Read the grid id = Dataset(grid_path, 'r') - lon = id.variables['lon_rho'][:,:] - lat = id.variables['lat_rho'][:,:] - mask_rho = id.variables['mask_rho'][:,:] - mask_zice = id.variables['mask_zice'][:,:] - zice = id.variables['zice'][:,:] + lon = id.variables['lon_rho'][:-15,:-2] + lat = id.variables['lat_rho'][:-15,:-2] + mask_rho = id.variables['mask_rho'][:-15,:-2] + mask_zice = id.variables['mask_zice'][:-15,:-2] + zice = id.variables['zice'][:-15,:-2] id.close() # Make sure longitude goes from -180 to 180, not 0 to 360 @@ -134,7 +134,7 @@ def massloss_map (grid_path, log_path, save=False, fig_name=None): # Determine bounds on colour scale max_val = amax(abs(error)) lev = linspace(-max_val, max_val, num=40) - lev = linspace(-200, 200, num=40) + lev = linspace(-100, 100, num=40) # Space ticks on colorbar 25% apart max_tick = floor(max_val/25)*25 diff --git a/max_vel.py b/max_vel.py index 0695d1e..3d19b36 100644 --- a/max_vel.py +++ b/max_vel.py @@ -17,8 +17,8 @@ def max_vel (file_path): for l in range(size(time)): print 'Processing timestep ' + str(l+1) + ' of ' + str(size(time)) - u = file.variables['u'][l,:,:,:] - v = file.variables['v'][l,:,:,:] + u = file.variables['u'][l,:,:-15,:-2] + v = file.variables['v'][l,:,:-15,:-3] # Find the maximum |u| and |v| at this timestep max_u.append(amax(abs(u))) max_v.append(amax(abs(v))) diff --git a/nsidc_aice_monthly.py b/nsidc_aice_monthly.py index 40c4bba..95adfc8 100644 --- a/nsidc_aice_monthly.py +++ b/nsidc_aice_monthly.py @@ -30,8 +30,8 @@ def nsidc_aice_monthly (cice_file, month, save=False, fig_name=None): # Read CICE grid and time values id = Dataset(cice_file, 'r') - cice_lon = id.variables['TLON'][:,:] - cice_lat = id.variables['TLAT'][:,:] + cice_lon = id.variables['TLON'][:-15,:] + cice_lat = id.variables['TLAT'][:-15,:] time_id = id.variables['time'] # Get the year, month, and day (all 1-based) for each output step # These are 5-day averages marked with the last day's date. @@ -109,12 +109,12 @@ def nsidc_aice_monthly (cice_file, month, save=False, fig_name=None): return # Start accumulating data weighted by days - cice_data[:,:] += id.variables['aice'][start_t,:,:]*start_days + cice_data[:,:] += id.variables['aice'][start_t,:-15,:]*start_days num_days += start_days # Beween start_t and end_t, we want all the days for t in range(start_t+1, end_t): - cice_data[:,:] += id.variables['aice'][t,:,:]*5 + cice_data[:,:] += id.variables['aice'][t,:-15,:]*5 num_days += 5 # Figure out how many of the 5 days averaged in end_t are actually within @@ -139,7 +139,7 @@ def nsidc_aice_monthly (cice_file, month, save=False, fig_name=None): print 'Error: ending index is month ' + str(cice_time[end_t].month) + ', day ' + str(cice_time[end_t].day) return - cice_data[:,:] += id.variables['aice'][end_t,:,:]*end_days + cice_data[:,:] += id.variables['aice'][end_t,:-15,:]*end_days num_days += end_days # Check that we got the correct number of days diff --git a/nsidc_aice_seasonal.py b/nsidc_aice_seasonal.py index 3811f79..efbed78 100644 --- a/nsidc_aice_seasonal.py +++ b/nsidc_aice_seasonal.py @@ -36,8 +36,8 @@ def nsidc_aice_seasonal (cice_file, save=False, fig_name=None): # Read CICE grid and time values id = Dataset(cice_file, 'r') - cice_lon = id.variables['TLON'][:,:] - cice_lat = id.variables['TLAT'][:,:] + cice_lon = id.variables['TLON'][:-15,:] + cice_lat = id.variables['TLAT'][:-15,:] time_id = id.variables['time'] # Get the year, month, and day (all 1-based) for each output step # These are 5-day averages marked with the last day's date. @@ -143,12 +143,12 @@ def nsidc_aice_seasonal (cice_file, save=False, fig_name=None): return # Start accumulating data weighted by days - cice_data[season,:,:] += id.variables['aice'][start_t_season,:,:]*start_days + cice_data[season,:,:] += id.variables['aice'][start_t_season,:-15,:]*start_days season_days += start_days # Beween start_t_season and end_t_season, we want all the days for t in range(start_t_season+1, end_t_season): - cice_data[season,:,:] += id.variables['aice'][t,:,:]*5 + cice_data[season,:,:] += id.variables['aice'][t,:-15,:]*5 season_days += 5 # Figure out how many of the 5 days averaged in end_t_season are @@ -172,7 +172,7 @@ def nsidc_aice_seasonal (cice_file, save=False, fig_name=None): print 'Error for season ' + season_names[season] + ': ending index is month ' + str(cice_time[end_t_season].month) + ', day ' + str(cice_time[end_t_season].day) return - cice_data[season,:,:] += id.variables['aice'][end_t_season,:,:]*end_days + cice_data[season,:,:] += id.variables['aice'][end_t_season,:-15,:]*end_days season_days += end_days # Check that we got the correct number of days @@ -253,23 +253,23 @@ def nsidc_aice_seasonal (cice_file, save=False, fig_name=None): lev = linspace(0, 1, num=50) # Plot - fig = figure(figsize=(9,18.5)) + fig = figure(figsize=(20,9)) # Loop over seasons for season in range(4): # NSIDC - ax = fig.add_subplot(4, 2, 2*season+1, aspect='equal') + ax = fig.add_subplot(2, 4, season+1, aspect='equal') contourf(nsidc_x, nsidc_y, nsidc_data[season,:,:], lev) if season == 0: - title('NSIDC', fontsize=24) + text(-39, 0, 'NSIDC', fontsize=24, ha='right') + title(season_names[season], fontsize=24) xlim([bdry1, bdry2]) ylim([bdry3, bdry4]) axis('off') # CICE - ax = fig.add_subplot(4, 2, 2*season+2, aspect='equal') + ax = fig.add_subplot(2, 4, season+5, aspect='equal') img = contourf(cice_x, cice_y, cice_data[season,:,:], lev) if season == 0: - title('CICE', fontsize=24) - text(35, 0, season_names[season], fontsize=24) + text(-39, 0, 'CICE', fontsize=24, ha='right') xlim([bdry1, bdry2]) ylim([bdry3, bdry4]) axis('off') diff --git a/other_models_data b/other_models_data new file mode 100644 index 0000000..9e2a262 --- /dev/null +++ b/other_models_data @@ -0,0 +1,43 @@ +FESOM (Timmermann 2012): +Name Mass loss Melt rate +Ronne-Filchner 138 0.35 +Brunt & Riiser-Larsen 65 0.94 +Fimbul & Jelbart 130 2.8 +Larsen C 48 1.0 +George VI 86 3.6 +Abbot 59 2.1 +Pine Island 13 3.1 +Getz 164 5.4 +Ross 260 0.6 +Amery 174 2.9 +Total 1600 + +BRIOS (Hellmer 2004): +Name Mass loss Melt rate +Abbot 18.16 0.55 +Amery 17.65 0.35 +Brunt & Riiser-Larsen 165.9 2.38 +Ronne-Filchner 119.7 0.32 +Fimbul 243.1 4.91 +George VI 22.48 0.43 +Getz 53.64 1.95 +Larsen (?) 38.13 0.63 +Ross 180.2 0.49 +Shackleton 47.68 1.04 +Total 906.6 + +COCO (Kusahara 2013): +Name Mass loss +Larsen (all) 21.6 +Ronne-Filchner 160.5 +Brunt to Prince Harald 107.6 +Amery 70.3 +West and Shackleton 27.1 +Totten to Merz 43.5 +McMurdo Sound 12.8 +Ross ` 110.7 +Sulzberger-Nickerson 4.8 +Getz to PIG 58.7 +Abbot 68.6 +Wilkins-GVI-Stange 83.3 +Total 769.6 \ No newline at end of file diff --git a/overturning_plot.py b/overturning_plot.py index 3840fe5..9ab64ea 100644 --- a/overturning_plot.py +++ b/overturning_plot.py @@ -18,12 +18,12 @@ def overturning_plot (file_path, fig_name): # Read v and grid variables id = Dataset(file_path, 'r') - v = id.variables['v'][-1,:,:,:] # -1 means last timestep - h = id.variables['h'][:,:] - zice = id.variables['zice'][:,:] - zeta = id.variables['zeta'][-1,:,:] - lon_v = id.variables['lon_v'][:,:] - lat_v = id.variables['lat_v'][:,:] + v = id.variables['v'][-1,:,:-15,:-3] # -1 means last timestep + h = id.variables['h'][:-15,:-3] + zice = id.variables['zice'][:-15,:-3] + zeta = id.variables['zeta'][-1,:-15,:-3] + lon_v = id.variables['lon_v'][:-15,:-3] + lat_v = id.variables['lat_v'][:-15,:-3] id.close() # Interpolate h, zice, and zeta to the v-grid @@ -60,7 +60,8 @@ def overturning_plot (file_path, fig_name): ylabel('Depth (m)') title('Meridional Overturning Streamfunction (Sv)') - savefig(fig_name) + #savefig(fig_name) + show() # Command-line interface diff --git a/sose_roms_seasonal.py b/sose_roms_seasonal.py index e9a2568..02a19a0 100644 --- a/sose_roms_seasonal.py +++ b/sose_roms_seasonal.py @@ -70,10 +70,10 @@ def sose_roms_seasonal (file_path, var_name, lon0, depth_bdry, save=False, fig_n # Read grid and time values id = Dataset(file_path, 'r') - h = id.variables['h'][:,:] - zice = id.variables['zice'][:,:] - lon_roms_2d = id.variables['lon_rho'][:,:] - lat_roms_2d = id.variables['lat_rho'][:,:] + h = id.variables['h'][:-15,:] + zice = id.variables['zice'][:-15,:] + lon_roms_2d = id.variables['lon_rho'][:-15,:] + lat_roms_2d = id.variables['lat_rho'][:-15,:] time_id = id.variables['ocean_time'] # Get the year, month, and day (all 1-based) for each output step # These are 5-day averages marked with the middle day's date @@ -186,12 +186,12 @@ def sose_roms_seasonal (file_path, var_name, lon0, depth_bdry, save=False, fig_n return # Start accumulating data weighted by days - var_3d_roms[season,:,:,:] += id.variables[var_name][start_t_season,:,:,:]*start_days + var_3d_roms[season,:,:,:] += id.variables[var_name][start_t_season,:,:-15,:]*start_days season_days += start_days # Between start_t_season and end_t_season, we want all the days for t in range(start_t_season+1, end_t_season): - var_3d_roms[season,:,:,:] += id.variables[var_name][t,:,:,:]*5 + var_3d_roms[season,:,:,:] += id.variables[var_name][t,:,:-15,:]*5 season_days += 5 # Figure out how many of the 5 days averaged in end_t_season are @@ -215,7 +215,7 @@ def sose_roms_seasonal (file_path, var_name, lon0, depth_bdry, save=False, fig_n print 'Error for season ' + season_title[season] + ': ending index is month ' + str(time[end_t_season].month) + ', day ' + str(time[end_t_season].day) return - var_3d_roms[season,:,:,:] += id.variables[var_name][end_t_season,:,:,:]*end_days + var_3d_roms[season,:,:,:] += id.variables[var_name][end_t_season,:,:-15,:]*end_days season_days += end_days # Check that we got the correct number of days diff --git a/spinup_plots.py b/spinup_plots.py index 7968b6c..2689bc1 100644 --- a/spinup_plots.py +++ b/spinup_plots.py @@ -35,11 +35,11 @@ def calc_grid (file_path): # 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] id.close() # Calculate water column thickness @@ -73,7 +73,7 @@ def get_rho (file_path, t): id = Dataset(file_path, 'r') # Read density anomalies, add rho0 to get absolute density # Convert to float128 to prevent overflow later - rho = ma.asarray(id.variables['rho'][t,:,:,:], dtype=float128) + rho0 + rho = ma.asarray(id.variables['rho'][t,:,:-15,:-3], dtype=float128) + rho0 id.close() return rho @@ -95,7 +95,7 @@ def calc_ohc (file_path, dV, rho, t): # Read temperature, converting to float128 to prevent overflow during # integration id = Dataset(file_path, 'r') - temp = ma.asarray(id.variables['temp'][t,:,:,:], dtype=float128) + temp = ma.asarray(id.variables['temp'][t,:,:-15,:-3], dtype=float128) # Convert from Celsius to Kelvin temp = temp + celsius2kelvin id.close() @@ -117,7 +117,7 @@ def calc_totalsalt (file_path, dV, rho, t): # Read salinity, converting to float128 to prevent overflow during # integration id = Dataset(file_path, 'r') - salt = ma.asarray(id.variables['salt'][t,:,:,:], dtype=float128) + salt = ma.asarray(id.variables['salt'][t,:,:-15,:-3], dtype=float128) id.close() # Integrate 1e-3*salt*rho over volume to get total mass of salt @@ -141,7 +141,7 @@ def calc_massloss (file_path, dA, t): # Read ice shelf melt rate, converting to float128 to prevent overflow # during integration id = Dataset(file_path, 'r') - ismr = ma.asarray(id.variables['m'][t,:,:], dtype=float128) + ismr = ma.asarray(id.variables['m'][t,:-15,:-3], dtype=float128) # Convert from m/s to m/y ismr = ismr*365.25*24*60*60 id.close() @@ -166,8 +166,8 @@ def calc_tke (file_path, dV, rho, t): # Read u and v, converting to float 128 to prevent overflow during # integration id = Dataset(file_path, 'r') - u = ma.asarray(id.variables['u'][t,:,:,:], dtype=float128) - v = ma.asarray(id.variables['v'][t,:,:,:], dtype=float128) + u = ma.asarray(id.variables['u'][t,:,:-15,:], dtype=float128) + v = ma.asarray(id.variables['v'][t,:,:-15,:], dtype=float128) id.close() # Interpolate u onto the rho-grid @@ -175,12 +175,16 @@ def calc_tke (file_path, dV, rho, t): middle_u = 0.5*(u[:,:,0:-1] + u[:,:,1:]) e_bdry_u = w_bdry_u[:,:] u_rho = ma.concatenate((w_bdry_u[:,:,None], middle_u, e_bdry_u[:,:,None]), axis=2) + # Throw away periodic boundary overlap + u_rho = u_rho[:,:,:-3] # Interpolate v onto the rho-grid s_bdry_v = v[:,0,:] middle_v = 0.5*(v[:,0:-1,:] + v[:,1:,:]) n_bdry_v = v[:,-1,:] v_rho = ma.concatenate((s_bdry_v[:,None,:], middle_v, n_bdry_v[:,None,:]), axis=1) + # Throw away periodic boundary overlap + v_rho = v_rho[:,:,:-3] # Integrate 0.5*rho*(u^2 + v^2) over volume to get TKE tke = sum(0.5*rho*(u_rho**2 + v_rho**2)*dV) @@ -211,7 +215,7 @@ def calc_drakepsgtrans (file_path, dy_wct, t): # Read ubar, converting to float128 to prevent overflow during integration id = Dataset(file_path, 'r') - ubar = ma.asarray(id.variables['ubar'][t,:,:], dtype=float128) + ubar = ma.asarray(id.variables['ubar'][t,:-15,:], dtype=float128) id.close() # Interpolate ubar onto the rho-grid @@ -219,6 +223,8 @@ def calc_drakepsgtrans (file_path, dy_wct, t): 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 periodic boundary overlap + ubar = ubar[:,:-3] # Trim arrays to these bounds ubar_rho_DP = ubar_rho[j_min:j_max,i_DP] @@ -241,7 +247,7 @@ def calc_totalice (cice_path, dA, t): id = Dataset(cice_path, 'r') # Read sea ice area fraction at each grid cell - aice = ma.asarray(id.variables['aice'][t,:,:], dtype=float128) + aice = ma.asarray(id.variables['aice'][t,:-15,:-3], dtype=float128) id.close() # Remove masks on aice and dA, and fill aice with zeros on land mask @@ -270,8 +276,8 @@ def calc_totalfwflux (cice_path, dA, t): id = Dataset(cice_path, 'r') # Read freshwater and salt fluxes at each grid cell - fresh = ma.asarray(id.variables['fresh_ai'][t,:,:], dtype=float128)/100./24./60./60 - fsalt = ma.asarray(id.variables['fsalt_ai'][t,:,:], dtype=float128)/1000. + fresh = ma.asarray(id.variables['fresh_ai'][t,:-15,:-3], dtype=float128)/100./24./60./60 + fsalt = ma.asarray(id.variables['fsalt_ai'][t,:-15,:-3], dtype=float128)/1000. fwflux = fresh - fsalt id.close() @@ -298,8 +304,8 @@ def calc_totalfwflux (cice_path, dA, t): def calc_bwtemp (file_path, dA, t): id = Dataset(file_path, 'r') - temp = ma.asarray(id.variables['temp'][t,0,:,:], dtype=float128) - zice = id.variables['zice'][:,:] + temp = ma.asarray(id.variables['temp'][t,0,:-15,:-3], dtype=float128) + zice = id.variables['zice'][:-15,:-3] id.close() cavity_temp = ma.masked_where(zice==0, temp) diff --git a/uv_vectorplot.py b/uv_vectorplot.py index 1febf17..50fbc9c 100644 --- a/uv_vectorplot.py +++ b/uv_vectorplot.py @@ -24,20 +24,20 @@ def uv_vectorplot (file_path, tstep, depth_key, save=False, fig_name=None): # Read grid and velocity data 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] if depth_key == 1: # Surface u and v - u = id.variables['u'][tstep-1,-1,:,:] - v = id.variables['v'][tstep-1,-1,:,:] + u = id.variables['u'][tstep-1,-1,:-15,:] + v = id.variables['v'][tstep-1,-1,:-15,:] elif depth_key == 2: # Bottom u and v - u = id.variables['u'][tstep-1,0,:,:] - v = id.variables['v'][tstep-1,0,:,:] + u = id.variables['u'][tstep-1,0,:-15,:] + v = id.variables['v'][tstep-1,0,:-15,:] elif depth_key == 3: # Vertically averaged u and v - u = id.variables['ubar'][tstep-1,:,:] - v = id.variables['vbar'][tstep-1,:,:] + u = id.variables['ubar'][tstep-1,:-15,:] + v = id.variables['vbar'][tstep-1,:-15,:] id.close() # Interpolate u to the rho-grid @@ -50,6 +50,9 @@ def uv_vectorplot (file_path, tstep, depth_key, save=False, fig_name=None): middle_v = 0.5*(v[0:-1,:] + v[1:,:]) n_bdry_v = v[-1,:] v_rho = ma.concatenate((s_bdry_v[None,:], middle_v, n_bdry_v[None,:]), axis=0) + # Throw away the overlapping periodic boundary + u_rho = u_rho[:,:-2] + v_rho = v_rho[:,:-2] # Calculate speed for the background filled contour plot speed = sqrt(u_rho**2 + v_rho**2) @@ -69,7 +72,7 @@ def uv_vectorplot (file_path, tstep, depth_key, save=False, fig_name=None): # Average X, Y, dX_dt, and dY_dt over block x block intervals # Calculate number of blocks size0 = int(ceil(size(X,0)/float(block))) - size1 = int(ceil(size(X,1)/float(block))) + size1 = int(ceil((size(X,1)-1)/float(block))) # Set up arrays for averaged fields X_block = ma.empty([size0, size1]) Y_block = ma.empty([size0, size1]) diff --git a/zice.py b/zice.py index fbb7458..af0a1ad 100644 --- a/zice.py +++ b/zice.py @@ -12,9 +12,9 @@ def zice (file_path, fig_name): # Read zice and masks file = Dataset(file_path, 'r') - zice = file.variables['zice'][:,:] - mask_rho = file.variables['mask_rho'][:,:] - mask_zice = file.variables['mask_zice'][:,:] + zice = file.variables['zice'][:-15,:-3] + mask_rho = file.variables['mask_rho'][:-15,:-3] + mask_zice = file.variables['mask_zice'][:-15,:-3] file.close() # Mask out ocean and land @@ -33,7 +33,8 @@ def zice (file_path, fig_name): xticks([814,815], (r'0$^{\circ}$',''), fontsize=24) yticks([34, 324], (r'75$^{\circ}$S', r'50$^{\circ}$S'), fontsize=24) - savefig(fig_name, transparent=True) + #savefig(fig_name, transparent=True) + show() # Command-line interface diff --git a/zice_circumpolar.py b/zice_circumpolar.py index a091e9c..71baf44 100644 --- a/zice_circumpolar.py +++ b/zice_circumpolar.py @@ -24,10 +24,10 @@ def zice_circumpolar (grid_path, fig_name): # Read data id = Dataset(grid_path, 'r') - zice = -1*id.variables['zice'][:,:] - lon = id.variables['lon_rho'][:,:] - lat = id.variables['lat_rho'][:,:] - mask_rho = id.variables['mask_rho'][:,:] + zice = -1*id.variables['zice'][:-15,:-2] + lon = id.variables['lon_rho'][:-15,:-2] + lat = id.variables['lat_rho'][:-15,:-2] + mask_rho = id.variables['mask_rho'][:-15,:-2] id.close() # Mask the open ocean and land @@ -68,8 +68,8 @@ def zice_circumpolar (grid_path, fig_name): title('Ice shelf draft (m)', fontsize=30) axis('off') - #show() - fig.savefig(fig_name) + show() + #fig.savefig(fig_name) # Command-line interface diff --git a/zonal_plot.py b/zonal_plot.py index b487d55..054c7af 100644 --- a/zonal_plot.py +++ b/zonal_plot.py @@ -33,9 +33,9 @@ def zonal_plot (file_path, var_name, tstep, lon_key, lon0, lon_bounds, depth_min # Read the variable id = Dataset(file_path, 'r') - data_3d = id.variables[var_name][tstep-1,:,:,:] + data_3d = id.variables[var_name][tstep-1,:,:-15,:] # Also read sea surface height - zeta = id.variables['zeta'][tstep-1,:,:] + zeta = id.variables['zeta'][tstep-1,:-15,:] if var_name == 'salt': units = 'psu' else: @@ -62,8 +62,8 @@ def zonal_plot (file_path, var_name, tstep, lon_key, lon0, lon_bounds, depth_min return # Read grid variables - h = id.variables['h'][:,:] - zice = id.variables['zice'][:,:] + h = id.variables['h'][:-15,:] + zice = id.variables['zice'][:-15,:] # h, zice, and zeta are on the rho-grid; interpolate if necessary if grid_name == 'u': h = 0.5*(h[:,0:-1] + h[:,1:]) @@ -74,10 +74,26 @@ def zonal_plot (file_path, var_name, tstep, lon_key, lon0, lon_bounds, depth_min zice = 0.5*(zice[0:-1,:] + zice[1:,:]) zeta = 0.5*(zeta[0:-1,:] + zeta[1:,:]) # Read the correct lat and lon for this grid (determined previously) - lon_2d = id.variables[lon_name][:,:] - lat_2d = id.variables[lat_name][:,:] + lon_2d = id.variables[lon_name][:-15,:] + lat_2d = id.variables[lat_name][:-15,:] id.close() + # Throw away periodic boundary overlap + if grid_name == 'u': + data_3d = data_3d[:,:,:-1] + zeta = zeta[:,:-1] + h = h[:,:-1] + zice = zice[:,:-1] + lon_2d = lon_2d[:,:-1] + lat_2d = lat_2d[:,:-1] + else: + data_3d = data_3d[:,:,:-2] + zeta = zeta[:,:-2] + h = h[:,:-2] + zice = zice[:,:-2] + lon_2d = lon_2d[:,:-2] + lat_2d = lat_2d[:,:-2] + # 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, zeta)