diff --git a/file_guide.txt b/file_guide.txt index 1ed36c0..ef9d923 100644 --- a/file_guide.txt +++ b/file_guide.txt @@ -1302,23 +1302,23 @@ mip_zonal_cavity_ts.py: Create one 2x2 plot for each major ice shelf: zonal file and the FESOM mesh directory. It will output one png for each ice shelf. -mip_mld_seasonal.py: Make a 4x2 plot showing seasonally averaged mixed layer - depth (defined as in Sallee et al 2013: depth at which - potential density is 0.03 kg/m^3 higher than at the - surface) comparing MetROMS (top) and FESOM (bottom). - To run: First clone my "fesomtools" GitHub repository and - replace '/short/y99/kaa561/fesomtools' (near the - top of the file) with the path to the cloned - repository on your system. Next, make seasonal - climatologies of 3D temperature and salnity for - both MetROMS and FESOM during the entire - simulation, using seasonal_climatology_roms.py (in - roms_tools) and seasonal_climatology.py (in - fesomtools). Once the files are ready, open python - or ipython and type "run mip_mld_seasonal.py". The - script will prompt you for paths to these files as - well as the ROMS grid file and the FESOM mesh - directory. It will output a png file. +mip_mld.py: Make a 4x2 plot showing summer (DJF) and winter (JJA) mixed layer + depth (defined as in Sallee et al 2013: depth at which potential + density is 0.03 kg/m^3 higher than at the surface) comparing + MetROMS (top) and FESOM (bottom) where each plot is shown zoomed + out to include the whole ACC, and zoomed in to the continental + shelf. + To run: First clone my "fesomtools" GitHub repository and replace + '/short/y99/kaa561/fesomtools' (near the top of the file) + with the path to the cloned repository on your system. + Next, make seasonal climatologies of 3D temperature and + salinity for both MetROMS and FESOM during the entire + simulation, using seasonal_climatology_roms.py (in + roms_tools) and seasonal_climatology.py (in fesomtools). + Once the files are ready, open python or ipython and type + "run mip_mld.py". The script will prompt you for paths to + these files as well as the ROMS grid file and the FESOM + mesh directory. It will output a png file. mip_dpt_calc.py: Calculate the mean Drake Passage transport over 2002-2016, as well as the linear trend, for MetROMS, low-res FESOM, and @@ -1330,11 +1330,10 @@ mip_dpt_calc.py: Calculate the mean Drake Passage transport over 2002-2016, as script will prompt you for the paths to the 3 log files. -mip_ts_distribution.py: Make a 2x1 plot of T/S distributions on the continental - shelf (defined by regions south of 60S with bathymetry - shallower than 60S, plus all ice shelf cavities) in - MetROMS (left) and FESOM (right). Include the surface - freezing point and density contours. +mip_ts_distribution.py: Make a 2x1 plot of T/S distributions south of 65S, + colour-coded based on depth, in MetROMS (left) and + FESOM (right). Include the surface freezing point and + density contours. To run: First clone my "fesomtools" GitHub repository and replace '/short/y99/kaa561/fesomtools' (near the top of the file) with the path to the @@ -1348,6 +1347,24 @@ mip_ts_distribution.py: Make a 2x1 plot of T/S distributions on the continental mesh directory, and the FESOM time-averaged file. It will output a png. +mip_drift_slices.py: Make a 4x2 plot showing zonal slices through 0E of + temperature (left) and salinity (right), in MetROMS (top) + and FESOM (bottom). The top 2x2 plots show the average + over the last year of simulation, while the bottom 2x2 + plots show anomalies with respect to the first year. + To run: First clone my "fesomtools" GitHub repository and + replace '/short/y99/kaa561/fesomtools' (near the + top of the file) with the path to the cloned + repository on your system. Next, make + time-averaged files of temperature and salinity + for each model over the first and last year (just + use ncra). Then open python or ipython and type + "run mip_drift_slices.py". The script will prompt + you for the paths to the ROMS grid file, the ROMS + time-averaged files, the FESOM mesh directory, and + the FESOM time-averaged files. It will output a + png. + ***UTILITY FUNCTIONS*** diff --git a/mip_drift_slices.py b/mip_drift_slices.py new file mode 100644 index 0000000..67ba713 --- /dev/null +++ b/mip_drift_slices.py @@ -0,0 +1,289 @@ +from netCDF4 import Dataset +from numpy import * +from matplotlib.pyplot import * +from matplotlib.patches import Polygon +from matplotlib.collections import PatchCollection +from calc_z import * +from interp_lon_roms import * +# Import FESOM scripts (have to modify path first) +import sys +sys.path.insert(0, '/short/y99/kaa561/fesomtools') +from fesom_grid import * +from fesom_sidegrid import * + +# Make a 4x2 plot showing zonal slices through 0E of temperature (left) and +# salinity (right), in MetROMS (top) and FESOM (bottom). The top 2x2 plots show +# the average over the last year of simulation, while the bottom 2x2 plots show +# anomalies with respect to the first year. +# Input: +# roms_grid = path to ROMS grid file +# roms_file_first, roms_file_last = paths to files containing ROMS temperature +# and salinity averaged over the first year +# and the last year respectively +# fesom_mesh_path = path to FESOM mesh directory +# fesom_file_first, fesom_file_last = paths to files containing FESOM +# temperature and salinity averaged over +# the first year and the last year +# respectively +def mip_drift_slices (roms_grid, roms_file_first, roms_file_last, fesom_mesh_path, fesom_file_first, fesom_file_last): + + # Longitude to plot (0E) + lon0 = 0 + # Bounds on plot + lat_min = -73 + lat_max = -30 + depth_min = -5300 + 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 + temp_anom = 2 + salt_min = 33.9 + salt_max = 34.9 + salt_anom = 0.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 ROMS' + # Read grid variables we need + id = Dataset(roms_grid, 'r') + roms_lon_2d = id.variables['lon_rho'][:,:] + roms_lat_2d = id.variables['lat_rho'][:,:] + roms_h = id.variables['h'][:,:] + roms_zice = id.variables['zice'][:,:] + id.close() + # Read temperature and salinity + id = Dataset(roms_file_first, 'r') + roms_temp_first_3d = id.variables['temp'][0,:,:,:] + roms_salt_first_3d = id.variables['salt'][0,:,:,:] + id.close() + id = Dataset(roms_file_last, 'r') + roms_temp_last_3d = id.variables['temp'][0,:,:,:] + roms_salt_last_3d = id.variables['salt'][0,:,:,:] + id.close() + # Get anomalies for last year + roms_temp_anom_3d = roms_temp_last_3d - roms_temp_first_3d + roms_salt_anom_3d = roms_salt_last_3d - roms_salt_first_3d + # Get a 3D array of z-coordinates; sc_r and Cs_r are unused in this script + roms_z_3d, sc_r, Cs_r = calc_z(roms_h, roms_zice, theta_s, theta_b, hc, N) + # Make sure we are in the range 0-360 + if lon0 < 0: + lon0 += 360 + # Interpolate to lon0 + roms_temp, roms_z, roms_lat = interp_lon_roms(roms_temp_last_3d, roms_z_3d, roms_lat_2d, roms_lon_2d, lon0) + roms_salt, roms_z, roms_lat = interp_lon_roms(roms_salt_last_3d, roms_z_3d, roms_lat_2d, roms_lon_2d, lon0) + roms_temp_anom, roms_z, roms_lat = interp_lon_roms(roms_temp_anom_3d, roms_z_3d, roms_lat_2d, roms_lon_2d, lon0) + roms_salt_anom, roms_z, roms_lat = interp_lon_roms(roms_salt_anom_3d, roms_z_3d, roms_lat_2d, roms_lon_2d, lon0) + # Switch back to range -180-180 + if lon0 > 180: + lon0 -= 360 + + print 'Processing FESOM' + # Build regular elements + elements = fesom_grid(fesom_mesh_path) + # Read temperature and salinity + id = Dataset(fesom_file_first, 'r') + fesom_temp_first_nodes = id.variables['temp'][0,:] + fesom_salt_first_nodes = id.variables['salt'][0,:] + id.close() + id = Dataset(fesom_file_last, 'r') + fesom_temp_last_nodes = id.variables['temp'][0,:] + fesom_salt_last_nodes = id.variables['salt'][0,:] + id.close() + # Get anomalies for last year + fesom_temp_anom_nodes = fesom_temp_last_nodes - fesom_temp_first_nodes + fesom_salt_anom_nodes = fesom_salt_last_nodes - fesom_salt_first_nodes + # Make SideElements + selements_temp = fesom_sidegrid(elements, fesom_temp_last_nodes, lon0, lat_max) + selements_salt = fesom_sidegrid(elements, fesom_salt_last_nodes, lon0, lat_max) + selements_temp_anom = fesom_sidegrid(elements, fesom_temp_anom_nodes, lon0, lat_max) + selements_salt_anom = fesom_sidegrid(elements, fesom_salt_anom_nodes, 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: + # Make patch + coord = transpose(vstack((selm.y, selm.z))) + patches.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_anom = [] + for selm in selements_temp_anom: + fesom_temp_anom.append(selm.var) + fesom_salt_anom = [] + for selm in selements_salt_anom: + fesom_salt_anom.append(selm.var) + + # 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+300, 0+1000, 1000) + depth_labels = [] + for val in depth_ticks: + depth_labels.append(str(int(round(-val)))) + + print 'Plotting' + fig = figure(figsize=(16,25)) + # Absolute values + gs1 = GridSpec(2,2) + gs1.update(left=0.1, right=0.95, bottom=0.57, top=0.925, wspace=0.08, hspace=0.2) + # MetROMS temperature + ax = subplot(gs1[0,0]) + pcolor(roms_lat, roms_z, roms_temp, vmin=temp_min, vmax=temp_max, cmap='jet') + title(r'MetROMS temperature ($^{\circ}$)', 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([]) + ax.set_yticks(depth_ticks) + ax.set_yticklabels(depth_labels, fontsize=16) + text(-38, 1000, 'Last year, ' + lon_string, fontsize=30) + # MetROMS salinity + ax = subplot(gs1[0,1]) + pcolor(roms_lat, roms_z, roms_salt, vmin=salt_min, vmax=salt_max, cmap='jet') + title('MetROMS salinity (psu)', fontsize=24) + xlim([lat_min, lat_max]) + ylim([depth_min, depth_max]) + ax.set_xticks(lat_ticks) + ax.set_xticklabels([]) + ax.set_yticks(depth_ticks) + ax.set_yticklabels([]) + # FESOM temperature + ax = subplot(gs1[1,0]) + img = PatchCollection(patches, cmap='jet') + img.set_array(array(fesom_temp)) + img.set_edgecolor('face') + img.set_clim(vmin=temp_min, vmax=temp_max) + ax.add_collection(img) + title(r'FESOM (high-res) temperature ($^{\circ}$C)', fontsize=24) + xlabel('Latitude', fontsize=18) + 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) + # Add a colorbar for temperature + cbaxes = fig.add_axes([0.17, 0.51, 0.3, 0.02]) + cbar = colorbar(img, orientation='horizontal', cax=cbaxes, extend='both', ticks=arange(temp_min, temp_max+2, 2)) + cbar.ax.tick_params(labelsize=16) + # FESOM salinity + ax = subplot(gs1[1,1]) + img = PatchCollection(patches, cmap='jet') + img.set_array(array(fesom_salt)) + img.set_edgecolor('face') + img.set_clim(vmin=salt_min, vmax=salt_max) + ax.add_collection(img) + title('FESOM (high-res) salinity (psu)', fontsize=24) + xlabel('Latitude', 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([]) + # Add a colorbar for salinity + cbaxes = fig.add_axes([0.6, 0.51, 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) + # Anomalies for last year + gs2 = GridSpec(2,2) + gs2.update(left=0.1, right=0.95, bottom=0.075, top=0.43, wspace=0.08, hspace=0.2) + # MetROMS temperature + ax = subplot(gs2[0,0]) + pcolor(roms_lat, roms_z, roms_temp_anom, vmin=-temp_anom, vmax=temp_anom, cmap='RdBu_r') + title(r'MetROMS temperature ($^{\circ}$)', 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([]) + ax.set_yticks(depth_ticks) + ax.set_yticklabels(depth_labels, fontsize=16) + text(-47, 1000, 'Last year minus first year, ' + lon_string, fontsize=30) + # MetROMS salinity + ax = subplot(gs2[0,1]) + pcolor(roms_lat, roms_z, roms_salt_anom, vmin=-salt_anom, vmax=salt_anom, cmap='RdBu_r') + title('MetROMS salinity (psu)', fontsize=24) + xlim([lat_min, lat_max]) + ylim([depth_min, depth_max]) + ax.set_xticks(lat_ticks) + ax.set_xticklabels([]) + ax.set_yticks(depth_ticks) + ax.set_yticklabels([]) + # FESOM temperature + ax = subplot(gs2[1,0]) + img = PatchCollection(patches, cmap='RdBu_r') + img.set_array(array(fesom_temp_anom)) + img.set_edgecolor('face') + img.set_clim(vmin=-temp_anom, vmax=temp_anom) + ax.add_collection(img) + title(r'FESOM (high-res) temperature ($^{\circ}$C)', fontsize=24) + xlabel('Latitude', fontsize=18) + 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) + # Add a colorbar for temperature + cbaxes = fig.add_axes([0.17, 0.02, 0.3, 0.02]) + cbar = colorbar(img, orientation='horizontal', cax=cbaxes, extend='both', ticks=arange(-temp_anom, temp_anom+1, 1)) + cbar.ax.tick_params(labelsize=16) + # FESOM salinity + ax = subplot(gs2[1,1]) + img = PatchCollection(patches, cmap='RdBu_r') + img.set_array(array(fesom_salt_anom)) + img.set_edgecolor('face') + img.set_clim(vmin=-salt_anom, vmax=salt_anom) + ax.add_collection(img) + title('FESOM (high-res) salinity (psu)', fontsize=24) + xlabel('Latitude', 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([]) + # Add a colorbar for salinity + cbaxes = fig.add_axes([0.6, 0.02, 0.3, 0.02]) + cbar = colorbar(img, orientation='horizontal', cax=cbaxes, extend='both', ticks=arange(-salt_anom, salt_anom+0.25, 0.25)) + cbar.ax.tick_params(labelsize=16) + #fig.show() + fig.savefig('ts_slices.png') + + +# Command-line interface +if __name__ == "__main__": + + roms_grid = raw_input("Path to ROMS grid file: ") + roms_file_first = raw_input("Path to ROMS temperature/salinity file averaged over first year: ") + roms_file_last = raw_input("Path to ROMS temperature/salinity file averaged over last year: ") + fesom_mesh_path = raw_input("Path to FESOM mesh directory: ") + fesom_file_first = raw_input("Path to FESOM temperature/salinity file averaged over first year: ") + fesom_file_last = raw_input("Path to FESOM temperature/salinity file averaged over last year: ") + mip_drift_slices(roms_grid, roms_file_first, roms_file_last, fesom_mesh_path, fesom_file_first, fesom_file_last) + + + + + diff --git a/mip_mld.py b/mip_mld.py new file mode 100644 index 0000000..79f79ff --- /dev/null +++ b/mip_mld.py @@ -0,0 +1,249 @@ +from netCDF4 import Dataset +from numpy import * +from matplotlib.collections import PatchCollection +from matplotlib.pyplot import * +from calc_z import * +# Import FESOM scripts (have to modify path first) +import sys +sys.path.insert(0, '/short/y99/kaa561/fesomtools') +from patches import * +# This will use the FESOM version of unesco.py for both MetROMS and FESOM, +# luckily it's identical +from unesco import * + +# Make a 4x2 plot showing summer (DJF) and winter (JJA) mixed layer depth +# (defined as in Sallee et al 2013: depth at which potential density is 0.03 +# kg/m^3 higher than at the surface) comparing MetROMS (top) and FESOM (bottom) +# where each plot is shown zoomed out to include the whole ACC, and zoomed in +# to the continental shelf. +# Input: +# 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): + + # Definition of mixed layer depth: where potential density exceeds + # surface density by this amount (kg/m^3) as in Sallee et al 2013 + density_anom = 0.03 + # Northern boundary for ACC plot: 30S + nbdry1 = -30 + 90 + # Northern boundary for continental shelf plot: 63S + nbdry2 = -64 + 90 + # Degrees to radians conversion factor + deg2rad = pi/180.0 + # FESOM parameters + circumpolar = True + mask_cavities = False + # ROMS parameters + theta_s = 7.0 + theta_b = 2.0 + hc = 250 + N = 31 + # Season names + season_names = ['DJF', 'MAM', 'JJA', 'SON'] + # Maximum for colour scale in each season + max_bound_summer = 150 + max_bound_winter = 600 + + print 'Processing MetROMS:' + print 'Reading grid' + id = Dataset(roms_grid, 'r') + roms_h = id.variables['h'][:,:] + roms_zice = id.variables['zice'][:,:] + roms_lon = id.variables['lon_rho'][:,:] + roms_lat = id.variables['lat_rho'][:,:] + id.close() + # Polar coordinates for plotting + roms_x = -(roms_lat+90)*cos(roms_lon*deg2rad+pi/2) + roms_y = (roms_lat+90)*sin(roms_lon*deg2rad+pi/2) + # Get a 3D array of z-coordinates; sc_r and Cs_r are unused in this script + roms_z, sc_r, Cs_r = calc_z(roms_h, roms_zice, theta_s, theta_b, hc, N) + # Make depth positive + roms_z = -1*roms_z + print 'Reading data' + id = Dataset(roms_seasonal_file, 'r') + roms_temp = id.variables['temp'][:,:,:,:] + roms_salt = id.variables['salt'][:,:,:,:] + id.close() + print 'Calculating density' + roms_density = unesco(roms_temp, roms_salt, zeros(shape(roms_temp))) + print 'Calculating mixed layer depth' + roms_mld = ma.empty([4, size(roms_lon,0), size(roms_lon,1)]) + # Awful triple loop here, can't find a cleaner way + for season in range(4): + print '...' + season_names[season] + for j in range(size(roms_lon,0)): + for i in range(size(roms_lon,1)): + # Get surface density + density_sfc = roms_density[season,-1,j,i] + # Get surface depth (only nonzero in ice shelf cavities) + depth_sfc = roms_z[-1,j,i] + if density_sfc is ma.masked: + # Land + roms_mld[season,j,i] = ma.masked + else: + # Loop downward + k = size(roms_density,1)-2 + while True: + if k < 0: + # Reached the bottom + roms_mld[season,j,i] = roms_z[0,j,i]-depth_sfc + break + if roms_density[season,k,j,i] >= density_sfc + density_anom: + # Reached the critical density anomaly + roms_mld[season,j,i] = roms_z[k,j,i]-depth_sfc + break + k -= 1 + + print 'Processing FESOM:' + print 'Building mesh' + elements, patches = make_patches(fesom_mesh_path, 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.close() + print 'Calculating density' + fesom_density_nodes = unesco(fesom_temp_nodes, fesom_salt_nodes, zeros(shape(fesom_temp_nodes))) + print 'Calculating mixed layer depth' + # Set up array for mixed layer depth at each element, at each season + fesom_mld = zeros([4, len(elements)]) + # Loop over seasons and elements to fill these in + for season in range(4): + print '...' + season_names[season] + mld_season = [] + for elm in elements: + # 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] + # 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[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[season,:] = array(mld_season) + + print 'Plotting' + fig = figure(figsize=(19,9)) + # MetROMS, summer, ACC + ax = fig.add_subplot(2, 4, 1, aspect='equal') + pcolor(roms_x, roms_y, roms_mld[0,:,:], vmin=0, vmax=max_bound_summer, cmap='jet') + text(-62, 0, 'MetROMS', fontsize=24, ha='right') + title(season_names[0], fontsize=24) + xlim([-nbdry1, nbdry1]) + ylim([-nbdry1, nbdry1]) + axis('off') + # MetROMS, summer, continental shelf + ax = fig.add_subplot(2, 4, 2, aspect='equal') + pcolor(roms_x, roms_y, roms_mld[0,:,:], vmin=0, vmax=max_bound_summer, cmap='jet') + title(season_names[0], fontsize=24) + xlim([-nbdry2, nbdry2]) + ylim([-nbdry2, nbdry2]) + axis('off') + # MetROMS, winter, ACC + ax = fig.add_subplot(2, 4, 3, aspect='equal') + pcolor(roms_x, roms_y, roms_mld[2,:,:], vmin=0, vmax=max_bound_winter, cmap='jet') + title(season_names[2], fontsize=24) + xlim([-nbdry1, nbdry1]) + ylim([-nbdry1, nbdry1]) + axis('off') + # MetROMS, winter, continental shelf + ax = fig.add_subplot(2, 4, 4, aspect='equal') + pcolor(roms_x, roms_y, roms_mld[2,:,:], vmin=0, vmax=max_bound_winter, cmap='jet') + title(season_names[2], fontsize=24) + xlim([-nbdry2, nbdry2]) + ylim([-nbdry2, nbdry2]) + axis('off') + # FESOM, summer, ACC + ax = fig.add_subplot(2, 4, 5, aspect='equal') + img = PatchCollection(patches, cmap='jet') + img.set_array(fesom_mld[0,:]) + img.set_clim(vmin=0, vmax=max_bound_summer) + img.set_edgecolor('face') + ax.add_collection(img) + xlim([-nbdry1, nbdry1]) + ylim([-nbdry1, nbdry1]) + axis('off') + text(-62, 0, 'FESOM', fontsize=24, ha='right') + text(-62, -10, '(high-res)', fontsize=24, ha='right') + # FESOM, summer, continental shelf + ax = fig.add_subplot(2, 4, 6, aspect='equal') + img = PatchCollection(patches, cmap='jet') + img.set_array(fesom_mld[0,:]) + img.set_clim(vmin=0, vmax=max_bound_summer) + img.set_edgecolor('face') + ax.add_collection(img) + xlim([-nbdry2, nbdry2]) + ylim([-nbdry2, nbdry2]) + axis('off') + # Add a horizontal colorbar for summer + cbaxes = fig.add_axes([0.2, 0.04, 0.2, 0.02]) + cbar = colorbar(img, orientation='horizontal', cax=cbaxes, extend='max', ticks=arange(0, max_bound_summer+50, 50)) + cbar.ax.tick_params(labelsize=20) + # FESOM, winter, ACC + ax = fig.add_subplot(2, 4, 7, aspect='equal') + img = PatchCollection(patches, cmap='jet') + img.set_array(fesom_mld[2,:]) + img.set_clim(vmin=0, vmax=max_bound_winter) + img.set_edgecolor('face') + ax.add_collection(img) + xlim([-nbdry1, nbdry1]) + ylim([-nbdry1, nbdry1]) + axis('off') + # FESOM, winter, continental shelf + ax = fig.add_subplot(2, 4, 8, aspect='equal') + img = PatchCollection(patches, cmap='jet') + img.set_array(fesom_mld[2,:]) + img.set_clim(vmin=0, vmax=max_bound_winter) + img.set_edgecolor('face') + ax.add_collection(img) + xlim([-nbdry2, nbdry2]) + ylim([-nbdry2, nbdry2]) + axis('off') + # Add a horizontal colorbar for winter + cbaxes = fig.add_axes([0.6, 0.04, 0.2, 0.02]) + cbar = colorbar(img, orientation='horizontal', cax=cbaxes, extend='max', ticks=arange(0, max_bound_winter+200, 200)) + cbar.ax.tick_params(labelsize=20) + # Add the main title + suptitle('Mixed layer depth (m), 2002-2016 average', fontsize=30) + # Decrease space between plots + subplots_adjust(wspace=0.025,hspace=0.025) + #fig.show() + fig.savefig('mld.png') + + +# Command-line interface +if __name__ == "__main__": + + 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) + + + + + + + + + diff --git a/mip_mld_seasonal.py b/mip_mld_seasonal.py deleted file mode 100644 index 644af08..0000000 --- a/mip_mld_seasonal.py +++ /dev/null @@ -1,202 +0,0 @@ -from netCDF4 import Dataset -from numpy import * -from matplotlib.collections import PatchCollection -from matplotlib.pyplot import * -from calc_z import * -# Import FESOM scripts (have to modify path first) -import sys -sys.path.insert(0, '/short/y99/kaa561/fesomtools') -from patches import * -# This will use the FESOM version of unesco.py for both MetROMS and FESOM, -# luckily it's identical -from unesco import * - -# Make a 4x2 plot showing seasonally averaged mixed layer depth (defined as in -# Sallee et al 2013: depth at which potential density is 0.03 kg/m^3 higher than -# at the surface) comparing MetROMS (top) and FESOM (bottom). -# Input: -# 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_seasonal (roms_grid, roms_seasonal_file, fesom_mesh_path, fesom_seasonal_file): - - # Definition of mixed layer depth: where potential density exceeds - # surface density by this amount (kg/m^3) as in Sallee et al 2013 - density_anom = 0.03 - # Northern boundary of plot 50S - nbdry = -50 + 90 - # Degrees to radians conversion factor - deg2rad = pi/180.0 - # FESOM parameters - circumpolar = True - mask_cavities = True - # ROMS parameters - theta_s = 7.0 - theta_b = 2.0 - hc = 250 - N = 31 - # Season names for plot titles - season_names = ['DJF', 'MAM', 'JJA', 'SON'] - # Maximum for colour scale - max_bound = 500 - - print 'Processing MetROMS:' - print 'Reading grid' - id = Dataset(roms_grid, 'r') - roms_h = id.variables['h'][:,:] - roms_zice = id.variables['zice'][:,:] - roms_lon = id.variables['lon_rho'][:,:] - roms_lat = id.variables['lat_rho'][:,:] - id.close() - # Polar coordinates for plotting - roms_x = -(roms_lat+90)*cos(roms_lon*deg2rad+pi/2) - roms_y = (roms_lat+90)*sin(roms_lon*deg2rad+pi/2) - # Get a 3D array of z-coordinates; sc_r and Cs_r are unused in this script - roms_z, sc_r, Cs_r = calc_z(roms_h, roms_zice, theta_s, theta_b, hc, N) - # Make depth positive - roms_z = -1*roms_z - print 'Reading data' - id = Dataset(roms_seasonal_file, 'r') - roms_temp = id.variables['temp'][:,:,:,:] - roms_salt = id.variables['salt'][:,:,:,:] - id.close() - print 'Calculating density' - roms_density = unesco(roms_temp, roms_salt, zeros(shape(roms_temp))) - print 'Calculating mixed layer depth' - roms_mld = ma.empty([4, size(roms_lon,0), size(roms_lon,1)]) - # Awful triple loop here, can't find a cleaner way - for season in range(4): - print '...' + season_names[season] - for j in range(size(roms_lon,0)): - for i in range(size(roms_lon,1)): - # Get surface density - density_sfc = roms_density[season,-1,j,i] - if density_sfc is ma.masked: - # Land - roms_mld[season,j,i] = ma.masked - elif roms_zice[j,i] != 0: - # Ice shelf - roms_mld[season,j,i] = ma.masked - else: - # Loop downward - k = size(roms_density,1)-2 - while True: - if k < 0: - # Reached the bottom - roms_mld[season,j,i] = roms_z[0,j,i] - break - if roms_density[season,k,j,i] >= density_sfc + density_anom: - # Reached the critical density anomaly - roms_mld[season,j,i] = roms_z[k,j,i] - break - k -= 1 - - print 'Processing FESOM:' - print 'Building mesh' - elements, patches = make_patches(fesom_mesh_path, 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.close() - print 'Calculating density' - fesom_density_nodes = unesco(fesom_temp_nodes, fesom_salt_nodes, zeros(shape(fesom_temp_nodes))) - print 'Calculating mixed layer depth' - # Count the number of elements not in ice shelf cavities - num_elm = 0 - for elm in elements: - if not elm.cavity: - num_elm += 1 - # Set up array for mixed layer depth at each element, at each season - fesom_mld = zeros([4, num_elm]) - # Loop over seasons and elements to fill these in - for season in range(4): - print '...' + season_names[season] - mld_season = [] - for elm in elements: - if not elm.cavity: - # Get mixed layer depth at each node - mld_nodes = [] - # Make sure we exclude ice shelf cavity nodes from element mean - # (an Element can be a non-cavity element and still have up to - # 2 cavity nodes) - for i in range(3): - if not elm.cavity_nodes[i]: - node = elm.nodes[i] - density_sfc = fesom_density_nodes[season,node.id] - temp_depth = node.depth - curr_node = node.below - while True: - if curr_node is None: - # Reached the bottom - mld_nodes.append(temp_depth) - break - if fesom_density_nodes[season,curr_node.id] >= density_sfc + density_anom: - # Reached the critical density anomaly - mld_nodes.append(curr_node.depth) - break - temp_depth = curr_node.depth - curr_node = curr_node.below - # For this element, save the mean mixed layer depth across - # non-cavity nodes (up to 3) - mld_season.append(mean(array(mld_nodes))) - fesom_mld[season,:] = array(mld_season) - - print 'Plotting' - fig = figure(figsize=(20,9)) - # Loop over seasons - for season in range(4): - # MetROMS - ax = fig.add_subplot(2, 4, season+1, aspect='equal') - pcolor(roms_x, roms_y, roms_mld[season,:,:], vmin=0, vmax=max_bound, cmap='jet') - if season == 0: - text(-43, 0, 'MetROMS', fontsize=24, ha='right') - title(season_names[season], fontsize=24) - xlim([-nbdry, nbdry]) - ylim([-nbdry, nbdry]) - axis('off') - # FESOM - ax = fig.add_subplot(2, 4, season+5, aspect='equal') - img = PatchCollection(patches, cmap='jet') - img.set_array(fesom_mld[season,:]) - img.set_clim(vmin=0, vmax=max_bound) - img.set_edgecolor('face') - ax.add_collection(img) - xlim([-nbdry, nbdry]) - ylim([-nbdry, nbdry]) - axis('off') - if season == 0: - text(-43, 0, 'FESOM', fontsize=24, ha='right') - # Add a horizontal colorbar at the bottom - cbaxes = fig.add_axes([0.25, 0.04, 0.5, 0.02]) - cbar = colorbar(img, orientation='horizontal', cax=cbaxes, extend='max') - cbar.ax.tick_params(labelsize=16) - # Add the main title - suptitle('Mixed layer depth (m)', fontsize=30) - # Decrease space between plots - subplots_adjust(wspace=0.025,hspace=0.025) - #fig.show() - fig.savefig('mld_seasonal.png') - - -# Command-line interface -if __name__ == "__main__": - - 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_seasonal (roms_grid, roms_seasonal_file, fesom_mesh_path, fesom_seasonal_file) - - - - - - - - - diff --git a/mip_ts_distribution.py b/mip_ts_distribution.py index 1fda29c..e124d9d 100644 --- a/mip_ts_distribution.py +++ b/mip_ts_distribution.py @@ -1,18 +1,16 @@ from netCDF4 import Dataset from numpy import * from matplotlib.pyplot import * +from matplotlib.colors import * from cartesian_grid_3d import * # Import FESOM scripts (have to modify path first) import sys sys.path.insert(0, '/short/y99/kaa561/fesomtools') from fesom_grid import * -# This will use the FESOM version of unesco.py for both MetROMS and FESOM, -# luckily it's identical from unesco import * -# Make a 2x1 plot of T/S distributions on the continental shelf (defined by -# regions south of 60S with bathymetry shallower than 60S, plus all ice shelf -# cavities) in MetROMS (left) and FESOM (right). Include the surface freezing +# Make a 2x1 plot of T/S distributions south of 65S, colour-coded based on +# depth, in MetROMS (left) and FESOM (right). Include the surface freezing # point and density contours. # Input: # roms_grid = path to ROMS grid file @@ -23,21 +21,20 @@ # salinity, over the same period as roms_file def mip_ts_distribution (roms_grid, roms_file, fesom_mesh_path, fesom_file): - # Bounds on latitude and bathymetry for continental shelf - lat0 = -60 - h0 = 1500 + # Northern boundary of water masses to consider + nbdry = -65 # Number of temperature and salinity bins num_bins = 1000 # Bounds on temperature and salinity bins (pre-computed, change if needed) min_salt = 32.3 - max_salt = 35.2 + max_salt = 35.1 min_temp = -3.1 - max_temp = 1.8 + max_temp = 3.8 # Bounds to actually plot min_salt_plot = 33.25 max_salt_plot = 35.0 min_temp_plot = -3 - max_temp_plot = 1.5 + max_temp_plot = 3.8 # FESOM grid generation parameters circumpolar = False cross_180 = False @@ -55,10 +52,13 @@ def mip_ts_distribution (roms_grid, roms_file, fesom_mesh_path, fesom_file): # Repeat for salinity salt_bins = linspace(min_salt, max_salt, num=num_bins) salt_centres = 0.5*(salt_bins[:-1] + salt_bins[1:]) - # Set up 2D arrays of temperature bins x salinity bins to increment with - # volume of water masses + # 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)]) + # Also arrays to integrate volume + volume_roms = zeros([size(temp_centres), size(salt_centres)]) + volume_fesom = 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) @@ -69,7 +69,7 @@ def mip_ts_distribution (roms_grid, roms_file, fesom_mesh_path, fesom_file): # salinity bins density = unesco(temp_2d, salt_2d, zeros(shape(temp_centres)))-1000 # Density contours to plot - density_lev = arange(26.8, 28.2, 0.2) + density_lev = arange(26.6, 28.4, 0.2) print 'Processing ROMS' # Read ROMS grid variables we need @@ -82,9 +82,9 @@ def mip_ts_distribution (roms_grid, roms_file, fesom_mesh_path, fesom_file): num_lat = size(roms_lat, 0) num_lon = size(roms_lon, 1) # Get integrands on 3D grid - dx, dy, dz, z = cartesian_grid_3d(roms_lon, roms_lat, roms_h, roms_zice, theta_s, theta_b, hc, N) + 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 = dx*dy*dz + dV = roms_dx*roms_dy*roms_dz # Read ROMS output id = Dataset(roms_file, 'r') roms_temp = id.variables['temp'][0,:,:,:] @@ -97,16 +97,20 @@ def mip_ts_distribution (roms_grid, roms_file, fesom_mesh_path, fesom_file): if roms_temp[0,j,i] is ma.masked: continue # Check if we're in the region of interest - if (roms_lat[j,i] < lat0 and roms_h[j,i] < h0) or (roms_zice[j,i] != 0): + if roms_lat[j,i] < nbdry: # Loop downward for k in range(N): # Figure out which bins this falls into temp_index = nonzero(temp_bins > roms_temp[k,j,i])[0][0] - 1 salt_index = nonzero(salt_bins > roms_salt[k,j,i])[0][0] - 1 - # Increment bins with volume - ts_vals_roms[temp_index, salt_index] += dV[k,j,i] + # Integrate depth*dV in this bin + ts_vals_roms[temp_index, salt_index] += -roms_z[k,j,i]*dV[k,j,i] + volume_roms[temp_index, salt_index] += dV[k,j,i] # Mask bins with zero volume - ts_vals_roms = ma.masked_where(ts_vals_roms==0, ts_vals_roms) + ts_vals_roms = ma.masked_where(volume_roms==0, ts_vals_roms) + volume_roms = ma.masked_where(volume_roms==0, volume_roms) + # Convert depths from integrals to volume-averages + ts_vals_roms /= volume_roms print 'Processing FESOM' # Make FESOM grid elements @@ -118,10 +122,8 @@ def mip_ts_distribution (roms_grid, roms_file, fesom_mesh_path, fesom_file): id.close() # Loop over elements for elm in elements: - # Find bathymetry at each node - node_bathy = array([elm.nodes[0].find_bottom().depth, elm.nodes[1].find_bottom().depth, elm.nodes[2].find_bottom().depth]) # See if we're in the region of interest - if (all(elm.lat < lat0) and all(node_bathy < h0)) or (elm.cavity): + if all(elm.lat < nbdry): # Get area of 2D triangle area = elm.area() nodes = [elm.nodes[0], elm.nodes[1], elm.nodes[2]] @@ -130,10 +132,11 @@ def mip_ts_distribution (roms_grid, roms_file, fesom_mesh_path, fesom_file): if nodes[0].below is None or nodes[1].below is None or nodes[2].below is None: # We've reached the bottom break - # Calculate average temperature, salinity, and layer thickness - # over this 3D triangular prism + # Calculate average temperature, salinity, depth, and layer + # thickness over this 3D triangular prism temp_vals = [] salt_vals = [] + depth_vals = [] dz = [] for i in range(3): # Average temperature over 6 nodes @@ -142,33 +145,43 @@ def mip_ts_distribution (roms_grid, roms_file, fesom_mesh_path, fesom_file): # Average salinity over 6 nodes salt_vals.append(fesom_salt[nodes[i].id]) salt_vals.append(fesom_salt[nodes[i].below.id]) + # Average depth over 6 nodes + depth_vals.append(nodes[i].depth) + depth_vals.append(nodes[i].below.depth) # Average dz over 3 vertical edges dz.append(abs(nodes[i].depth - nodes[i].below.depth)) # Get ready for next repetition of loop nodes[i] = nodes[i].below temp_elm = mean(array(temp_vals)) salt_elm = mean(array(salt_vals)) + depth_elm = mean(array(depth_vals)) # Calculate volume of 3D triangular prism volume = area*mean(array(dz)) # Figure out which bins this falls into temp_index = nonzero(temp_bins > temp_elm)[0][0] - 1 salt_index = nonzero(salt_bins > salt_elm)[0][0] - 1 - # Increment bins with volume - ts_vals_fesom[temp_index, salt_index] += volume + # Integrate depth*volume in this bin + ts_vals_fesom[temp_index, salt_index] += depth_elm*volume + volume_fesom[temp_index, salt_index] += volume # Mask bins with zero volume - ts_vals_fesom = ma.masked_where(ts_vals_fesom==0, ts_vals_fesom) + ts_vals_fesom = ma.masked_where(volume_fesom==0, ts_vals_fesom) + volume_fesom = ma.masked_where(volume_fesom==0, volume_fesom) + # Convert depths from integrals to volume-averages + ts_vals_fesom /= volume_fesom - # Find bounds on log of volume - min_val = min(amin(log(ts_vals_roms)), amin(log(ts_vals_fesom))) - max_val = max(amax(log(ts_vals_roms)), amax(log(ts_vals_fesom))) + # Find the maximum depth for plotting + max_depth = max(amax(ts_vals_roms), amax(ts_vals_fesom)) + # 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, 0.5), (33.6, 0.75), (33.9, 1.0), (34.2, 1.0), (34.45, 1.3), (34.75, 1.3), (34.9, 0.5)] + 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)] print "Plotting" fig = figure(figsize=(20,12)) # ROMS ax = fig.add_subplot(1, 2, 1) - pcolor(salt_centres, temp_centres, log(ts_vals_roms), vmin=min_val, vmax=max_val, cmap='jet') + 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') # Add density contours @@ -183,7 +196,7 @@ def mip_ts_distribution (roms_grid, roms_file, fesom_mesh_path, fesom_file): title('MetROMS', fontsize=28) # FESOM ax = fig.add_subplot(1, 2, 2) - img = pcolor(salt_centres, temp_centres, log(ts_vals_fesom), vmin=min_val, vmax=max_val, cmap='jet') + img = pcolor(salt_centres, temp_centres, ts_vals_fesom, 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) @@ -195,10 +208,10 @@ def mip_ts_distribution (roms_grid, roms_file, fesom_mesh_path, fesom_file): title('FESOM (high-res)', fontsize=28) # Add a colourbar on the right cbaxes = fig.add_axes([0.93, 0.2, 0.02, 0.6]) - cbar = colorbar(img, cax=cbaxes) + cbar = colorbar(img, cax=cbaxes, ticks=[0,100,500,1000,2000,3000,4000]) cbar.ax.tick_params(labelsize=18) # Add the main title - suptitle('Continental shelf water masses (log of volume), 2002-2016 average', fontsize=30) + suptitle('Water masses south of 65$^{\circ}$S: depth (m), 2002-2016 average', fontsize=30) subplots_adjust(wspace=0.1) #fig.show() fig.savefig('ts_distribution.png') @@ -212,8 +225,9 @@ def mip_ts_distribution (roms_grid, roms_file, fesom_mesh_path, fesom_file): 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) - - + + +