diff --git a/add_tide_period.py b/add_tide_period.py index 7c47048..3c552b0 100644 --- a/add_tide_period.py +++ b/add_tide_period.py @@ -8,7 +8,7 @@ def add_tide_period (): # Tide file to add periods to # Order is q1, o1, p1, k1, n2, m2, s2, k2 - tide_file = '../ROMS-CICE-MCT/data/pot_tides.nc' + tide_file = '../metroms_iceshelf/data/pot_tides.nc' # Periods in seconds # Order is m2, s2, n2, k2, k1, o1, p1, q1, mf, mm period = array([44714.165191868, 43200.0012869521, 86164.0770050671, 92949.6357005365, 45570.0535117177, 86637.1997716528, 43082.0503185947, 96726.0857029666, 2380715.86358729, 1180295.54554976]) diff --git a/adv_polynyas.py b/adv_polynyas.py index ce4c123..128aeef 100644 --- a/adv_polynyas.py +++ b/adv_polynyas.py @@ -2,7 +2,7 @@ from numpy import * from matplotlib.pyplot import * -# Create a 1x2 plot of sea ice concentration along the coast of Wilkes Land +# Create a 1x2 plot of sea ice concentration near the Amery Ice Shelf # on 23 August (the sea ice area max), showing the difference in polynya size # between the U3_LIM and C4_LD simulations. def adv_polynyas (): diff --git a/adv_ross_tsplots.py b/adv_ross_tsplots.py index a50788c..217ba66 100644 --- a/adv_ross_tsplots.py +++ b/adv_ross_tsplots.py @@ -93,8 +93,8 @@ def adv_ross_tsplots (): # Main title suptitle(r'180$^{\circ}$E, 31 December (Ross Sea)', fontsize=30) - fig.show() - #fig.savefig('adv_ross_tsplots.png') + #fig.show() + fig.savefig('adv_ross_tsplots.png') # Command-line interface diff --git a/salinity_diff.py b/adv_u3_sss_anom.py similarity index 90% rename from salinity_diff.py rename to adv_u3_sss_anom.py index 1e608f8..792973e 100644 --- a/salinity_diff.py +++ b/adv_u3_sss_anom.py @@ -2,10 +2,15 @@ from numpy import * from matplotlib.pyplot import * -def salinity_diff (): +# Plot the sea surface salinity anomaly between the U3 and U3_LIM simulations +# on 23 August (sea ice area maximum). +def adv_u3_sss_anom (): + # Paths to simulation directories paths = ['/short/m68/kaa561/advection/u3_lim/', '/short/m68/kaa561/advection/u3/'] + # File for 23 August daily average file_tail = 'iceh.1992-08-23.nc' + # Bounds on colour scale max_anom = 0.2 tick_anom = 0.1 # Degrees to radians conversion factor @@ -75,8 +80,9 @@ def salinity_diff (): fig.savefig('sss_u3_anom.png') +# Command-line interface if __name__ == "__main__": - salinity_diff() + adv_u3_sss_anom() diff --git a/aice_hs_seasonal.py b/aice_hs_seasonal.py deleted file mode 100644 index ea9c486..0000000 --- a/aice_hs_seasonal.py +++ /dev/null @@ -1,116 +0,0 @@ -from numpy import * -from netCDF4 import Dataset -from matplotlib.pyplot import * -from seasonal_avg_cice import * - -def aice_hs_seasonal (cice_file, save=False, fig_name=None): - - season_names = ['DJF', 'MAM', 'JJA', 'SON'] - deg2rad = pi/180.0 - - id = Dataset(cice_file, 'r') - lon_tmp = id.variables['TLON'][:-15,:] - lat_tmp = id.variables['TLAT'][:-15,:] - num_lon = id.variables['TLON'].shape[1] - num_lat = id.variables['TLAT'].shape[0] - id.close() - lon = ma.empty([size(lon_tmp,0), size(lon_tmp,1)+1]) - lat = ma.empty([size(lat_tmp,0), size(lat_tmp,1)+1]) - lon[:,:-1] = lon_tmp - lon[:,-1] = lon_tmp[:,0] - lat[:,:-1] = lat_tmp - lat[:,-1] = lat_tmp[:,0] - - aice_tmp = seasonal_avg_cice(cice_file, 'aice', [num_lat, num_lon]) - hs_tmp = seasonal_avg_cice(cice_file, 'hs', [num_lat, num_lon]) - aice_tmp = aice_tmp[:,:-15,:] - hs_tmp = hs_tmp[:,:-15,:] - aice = ma.empty([size(aice_tmp,0), size(aice_tmp,1), size(aice_tmp,2)+1]) - aice[:,:,:-1] = aice_tmp - aice[:,:,-1] = aice_tmp[:,:,0] - hs = ma.empty([size(hs_tmp,0), size(hs_tmp,1), size(hs_tmp,2)+1]) - hs[:,:,:-1] = hs_tmp - hs[:,:,-1] = hs_tmp[:,:,0] - - x = -(lat+90)*cos(lon*deg2rad+pi/2) - y = (lat+90)*sin(lon*deg2rad+pi/2) - - bdry1 = -35 - bdry2 = 35 - bdry3 = -33 - bdry4 = 37 - - # Set consistent colour levels - lev1 = linspace(0, 1, num=50) - lev2 = linspace(0, 0.5, num=50) - - # Plot - fig = figure(figsize=(20,9)) - # Loop over seasons - for season in range(4): - ax = fig.add_subplot(2, 4, season+1, aspect='equal') - img = contourf(x, y, aice[season,:,:], lev1) - xlim([bdry1, bdry2]) - ylim([bdry3, bdry4]) - axis('off') - if season == 0: - text(-39, 0, 'aice (1)', fontsize=21, ha='right') - title(season_names[season], fontsize=24) - if season == 3: - cbaxes1 = fig.add_axes([0.92, 0.55, 0.01, 0.3]) - cbar1 = colorbar(img, ticks=arange(0,1+0.25,0.25), cax=cbaxes1) - cbar1.ax.tick_params(labelsize=16) - ax = fig.add_subplot(2, 4, season+5, aspect='equal') - img = contourf(x, y, hs[season,:,:], lev2, extend='both') - xlim([bdry1, bdry2]) - ylim([bdry3, bdry4]) - axis('off') - if season == 0: - text(-39, 0, 'hs (m)', fontsize=21, ha='right') - if season == 3: - cbaxes2 = fig.add_axes([0.92, 0.15, 0.01, 0.3]) - cbar2 = colorbar(img, ticks=arange(0,0.5+0.1,0.1), cax=cbaxes2) - cbar2.ax.tick_params(labelsize=16) - subplots_adjust(wspace=0.025,hspace=0.025) - - # Finished - if save: - fig.savefig(fig_name) - else: - fig.show() - - -# Command-line interface -if __name__ == "__main__": - - cice_file = raw_input("Path to CICE file, containing at least one complete Dec-Nov period: ") - action = raw_input("Save figure (s) or display on screen (d)? ") - if action == 's': - save = True - fig_name = raw_input("File name for figure: ") - elif action == 'd': - save = False - fig_name = None - aice_hs_seasonal(cice_file, save, fig_name) - - - - - - - - - - - - - - - - - - - - - - diff --git a/average_erainterim.py b/average_erainterim.py deleted file mode 100644 index 07923d4..0000000 --- a/average_erainterim.py +++ /dev/null @@ -1,62 +0,0 @@ -from os import system - -def average_erainterim (): - - head_in = '../ROMS-CICE-MCT/data/ERA_Interim/originals/FC_' - tail_in = '_unlim_orig.nc' - head_out = '../ROMS-CICE-MCT/data/ERA_Interim/monthly/FC_' - tail_out = '_monthly_orig.nc' - steps_per_day = 2 - - days_per_month = [31,28,31,30,31,30,31,31,30,31,30,31] - first_year = 1992 - last_year = 2005 - leap_years = range(1976,last_year+1,4) - - for year in range(first_year,last_year+1): - posn = 0 - for month in range(12): - print('Processing month '+str(month+1)+' of '+str(year)) - if month == 1 and year in leap_years: - days = 29 - else: - days = days_per_month[month] - num_steps = days*an_steps_per_day - month_str = str(month+1) - if month < 9: - month_str = str(0) + month_str - in_file = an_head_in + str(year) + an_tail_in - tmp_file = 'tmp' + month_str + '.nc' - system('ncra -d time,' + str(posn) + ',' + str(posn+num_steps-1) + ' ' + in_file + ' ' + tmp_file) - posn = posn + num_steps - print('Concatenating '+str(year)) - out_file = an_head_out + str(year) + an_tail_out - system('ncrcat tmp*.nc ' + out_file) - system('rm tmp*.nc') - - for year in range(first_year,last_year+1): - posn = 0 - for month in range(12): - print('Processing month '+str(month+1)+' of '+str(year)) - if month == 1 and year in leap_years: - days = 29 - else: - days = days_per_month[month] - num_steps = days*fc_steps_per_day - month_str = str(month+1) - if month < 9: - month_str = str(0) + month_str - in_file = fc_head_in + str(year) + fc_tail_in - tmp_file = 'tmp' + month_str + '.nc' - system('ncra -d time,' + str(posn) + ',' + str(posn+num_steps-1) + ' ' + in_file + ' ' + tmp_file) - posn = posn + num_steps - print('Concatenating '+str(year)) - out_file = fc_head_out + str(year) + fc_tail_out - system('ncrcat tmp*.nc ' + out_file) - system('rm tmp*.nc') - - -if __name__ == "__main__": - - average_erainterim() - diff --git a/cice_ini.py b/cice_ini.py index 2ac5f3c..3710baa 100644 --- a/cice_ini.py +++ b/cice_ini.py @@ -2,48 +2,72 @@ from numpy import * from scipy.interpolate import griddata +# Make a CICE initial conditions file (non-standard option ice_ic='mask') based +# on NSIDC satellite observations of sea ice concentration. Wherever +# concentration exceeds 15%, set it to 100%; everywhere else, set it to 0%. +# The CICE code will read this "mask" and set a constant sea ice thickness of +# 1 m and snow thickness of 0.2 m. +# Input: +# nsidc_file = path to NSIDC file: CDR monthly average of sea ice concentration +# roms_grid = path to ROMS grid file +# out_file = desired path to output file def cice_ini (nsidc_file, roms_grid, out_file): + # Read the NSIDC grid and data id = Dataset(nsidc_file, 'r') nsidc_lon = id.variables['longitude'][:,:] nsidc_lat = id.variables['latitude'][:,:] nsidc_aice = id.variables['seaice_conc_monthly_cdr'][0,:,:] id.close() + # Make sure longitude goes from 0 to 360 to fit with ROMS convention index = nsidc_lon < 0 nsidc_lon[index] = nsidc_lon[index] + 360 + # Make sure sea ice concentration values go from 0 to 1 index = nsidc_aice < 0 nsidc_aice[index] = 0.0 index = nsidc_aice > 1 nsidc_aice[index] = 1.0 + # Read the ROMS grid id = Dataset(roms_grid, 'r') roms_lon = id.variables['lon_rho'][:,:] roms_lat = id.variables['lat_rho'][:,:] id.close() + # Make an nx2 array of the NSIDC coordinates, flattened points = empty([size(nsidc_lon), 2]) points[:,0] = ravel(nsidc_lon) points[:,1] = ravel(nsidc_lat) + # Also flatten the NSIDC data values = ravel(nsidc_aice) + # Now make an mx2 array of the ROMS coordinates xi = empty([size(roms_lon), 2]) xi[:,0] = ravel(roms_lon) xi[:,1] = ravel(roms_lat) + # Linear interpolation to the ROMS coordinates result = griddata(points, values, xi, method='linear', fill_value=-999) + # Also do a nearest-neighbour interpolation result2 = griddata(points, values, xi, method='nearest') + # Replace missing values (from land mask) in linear interpolation with + # nearest-neighbour interpolation so there are no coastal artifacts index = result==-999 result[index] = result2[index] + # Reshape to the correct dimensions aice = reshape(result, shape(roms_lon)) + # Cut-off concentration of 15% index = aice > 0.15 aice[index] = 1.0 aice[~index] = 0.0 + # Remove the edges (CICE only wants interior points) aice = aice[1:-1,1:-1] num_lon = size(aice, 1) num_lat = size(aice, 0) + # Write output out_id = Dataset(out_file, 'w') out_id.createDimension('xi_rho', num_lon) out_id.createDimension('eta_rho', num_lat) @@ -58,6 +82,7 @@ def cice_ini (nsidc_file, roms_grid, out_file): out_id.close() +# Command-line interface if __name__ == "__main__": nsidc_file = '/short/m68/kaa561/nsidc_aice/seaice_conc_monthly_sh_f11_199201_v02r00.nc' diff --git a/colormaps.py b/colormaps.py deleted file mode 100644 index d1f74a4..0000000 --- a/colormaps.py +++ /dev/null @@ -1,1058 +0,0 @@ -# New matplotlib colormaps by Nathaniel J. Smith, Stefan van der Walt, -# and (in the case of viridis) Eric Firing. -# -# This file and the colormaps in it are released under the CC0 license / -# public domain dedication. We would appreciate credit if you use or -# redistribute these colormaps, but do not impose any legal restrictions. -# -# To the extent possible under law, the persons who associated CC0 with -# mpl-colormaps have waived all copyright and related or neighboring rights -# to mpl-colormaps. -# -# You should have received a copy of the CC0 legalcode along with this -# work. If not, see . - -__all__ = ['magma', 'inferno', 'plasma', 'viridis'] - -_magma_data = [[0.001462, 0.000466, 0.013866], - [0.002258, 0.001295, 0.018331], - [0.003279, 0.002305, 0.023708], - [0.004512, 0.003490, 0.029965], - [0.005950, 0.004843, 0.037130], - [0.007588, 0.006356, 0.044973], - [0.009426, 0.008022, 0.052844], - [0.011465, 0.009828, 0.060750], - [0.013708, 0.011771, 0.068667], - [0.016156, 0.013840, 0.076603], - [0.018815, 0.016026, 0.084584], - [0.021692, 0.018320, 0.092610], - [0.024792, 0.020715, 0.100676], - [0.028123, 0.023201, 0.108787], - [0.031696, 0.025765, 0.116965], - [0.035520, 0.028397, 0.125209], - [0.039608, 0.031090, 0.133515], - [0.043830, 0.033830, 0.141886], - [0.048062, 0.036607, 0.150327], - [0.052320, 0.039407, 0.158841], - [0.056615, 0.042160, 0.167446], - [0.060949, 0.044794, 0.176129], - [0.065330, 0.047318, 0.184892], - [0.069764, 0.049726, 0.193735], - [0.074257, 0.052017, 0.202660], - [0.078815, 0.054184, 0.211667], - [0.083446, 0.056225, 0.220755], - [0.088155, 0.058133, 0.229922], - [0.092949, 0.059904, 0.239164], - [0.097833, 0.061531, 0.248477], - [0.102815, 0.063010, 0.257854], - [0.107899, 0.064335, 0.267289], - [0.113094, 0.065492, 0.276784], - [0.118405, 0.066479, 0.286321], - [0.123833, 0.067295, 0.295879], - [0.129380, 0.067935, 0.305443], - [0.135053, 0.068391, 0.315000], - [0.140858, 0.068654, 0.324538], - [0.146785, 0.068738, 0.334011], - [0.152839, 0.068637, 0.343404], - [0.159018, 0.068354, 0.352688], - [0.165308, 0.067911, 0.361816], - [0.171713, 0.067305, 0.370771], - [0.178212, 0.066576, 0.379497], - [0.184801, 0.065732, 0.387973], - [0.191460, 0.064818, 0.396152], - [0.198177, 0.063862, 0.404009], - [0.204935, 0.062907, 0.411514], - [0.211718, 0.061992, 0.418647], - [0.218512, 0.061158, 0.425392], - [0.225302, 0.060445, 0.431742], - [0.232077, 0.059889, 0.437695], - [0.238826, 0.059517, 0.443256], - [0.245543, 0.059352, 0.448436], - [0.252220, 0.059415, 0.453248], - [0.258857, 0.059706, 0.457710], - [0.265447, 0.060237, 0.461840], - [0.271994, 0.060994, 0.465660], - [0.278493, 0.061978, 0.469190], - [0.284951, 0.063168, 0.472451], - [0.291366, 0.064553, 0.475462], - [0.297740, 0.066117, 0.478243], - [0.304081, 0.067835, 0.480812], - [0.310382, 0.069702, 0.483186], - [0.316654, 0.071690, 0.485380], - [0.322899, 0.073782, 0.487408], - [0.329114, 0.075972, 0.489287], - [0.335308, 0.078236, 0.491024], - [0.341482, 0.080564, 0.492631], - [0.347636, 0.082946, 0.494121], - [0.353773, 0.085373, 0.495501], - [0.359898, 0.087831, 0.496778], - [0.366012, 0.090314, 0.497960], - [0.372116, 0.092816, 0.499053], - [0.378211, 0.095332, 0.500067], - [0.384299, 0.097855, 0.501002], - [0.390384, 0.100379, 0.501864], - [0.396467, 0.102902, 0.502658], - [0.402548, 0.105420, 0.503386], - [0.408629, 0.107930, 0.504052], - [0.414709, 0.110431, 0.504662], - [0.420791, 0.112920, 0.505215], - [0.426877, 0.115395, 0.505714], - [0.432967, 0.117855, 0.506160], - [0.439062, 0.120298, 0.506555], - [0.445163, 0.122724, 0.506901], - [0.451271, 0.125132, 0.507198], - [0.457386, 0.127522, 0.507448], - [0.463508, 0.129893, 0.507652], - [0.469640, 0.132245, 0.507809], - [0.475780, 0.134577, 0.507921], - [0.481929, 0.136891, 0.507989], - [0.488088, 0.139186, 0.508011], - [0.494258, 0.141462, 0.507988], - [0.500438, 0.143719, 0.507920], - [0.506629, 0.145958, 0.507806], - [0.512831, 0.148179, 0.507648], - [0.519045, 0.150383, 0.507443], - [0.525270, 0.152569, 0.507192], - [0.531507, 0.154739, 0.506895], - [0.537755, 0.156894, 0.506551], - [0.544015, 0.159033, 0.506159], - [0.550287, 0.161158, 0.505719], - [0.556571, 0.163269, 0.505230], - [0.562866, 0.165368, 0.504692], - [0.569172, 0.167454, 0.504105], - [0.575490, 0.169530, 0.503466], - [0.581819, 0.171596, 0.502777], - [0.588158, 0.173652, 0.502035], - [0.594508, 0.175701, 0.501241], - [0.600868, 0.177743, 0.500394], - [0.607238, 0.179779, 0.499492], - [0.613617, 0.181811, 0.498536], - [0.620005, 0.183840, 0.497524], - [0.626401, 0.185867, 0.496456], - [0.632805, 0.187893, 0.495332], - [0.639216, 0.189921, 0.494150], - [0.645633, 0.191952, 0.492910], - [0.652056, 0.193986, 0.491611], - [0.658483, 0.196027, 0.490253], - [0.664915, 0.198075, 0.488836], - [0.671349, 0.200133, 0.487358], - [0.677786, 0.202203, 0.485819], - [0.684224, 0.204286, 0.484219], - [0.690661, 0.206384, 0.482558], - [0.697098, 0.208501, 0.480835], - [0.703532, 0.210638, 0.479049], - [0.709962, 0.212797, 0.477201], - [0.716387, 0.214982, 0.475290], - [0.722805, 0.217194, 0.473316], - [0.729216, 0.219437, 0.471279], - [0.735616, 0.221713, 0.469180], - [0.742004, 0.224025, 0.467018], - [0.748378, 0.226377, 0.464794], - [0.754737, 0.228772, 0.462509], - [0.761077, 0.231214, 0.460162], - [0.767398, 0.233705, 0.457755], - [0.773695, 0.236249, 0.455289], - [0.779968, 0.238851, 0.452765], - [0.786212, 0.241514, 0.450184], - [0.792427, 0.244242, 0.447543], - [0.798608, 0.247040, 0.444848], - [0.804752, 0.249911, 0.442102], - [0.810855, 0.252861, 0.439305], - [0.816914, 0.255895, 0.436461], - [0.822926, 0.259016, 0.433573], - [0.828886, 0.262229, 0.430644], - [0.834791, 0.265540, 0.427671], - [0.840636, 0.268953, 0.424666], - [0.846416, 0.272473, 0.421631], - [0.852126, 0.276106, 0.418573], - [0.857763, 0.279857, 0.415496], - [0.863320, 0.283729, 0.412403], - [0.868793, 0.287728, 0.409303], - [0.874176, 0.291859, 0.406205], - [0.879464, 0.296125, 0.403118], - [0.884651, 0.300530, 0.400047], - [0.889731, 0.305079, 0.397002], - [0.894700, 0.309773, 0.393995], - [0.899552, 0.314616, 0.391037], - [0.904281, 0.319610, 0.388137], - [0.908884, 0.324755, 0.385308], - [0.913354, 0.330052, 0.382563], - [0.917689, 0.335500, 0.379915], - [0.921884, 0.341098, 0.377376], - [0.925937, 0.346844, 0.374959], - [0.929845, 0.352734, 0.372677], - [0.933606, 0.358764, 0.370541], - [0.937221, 0.364929, 0.368567], - [0.940687, 0.371224, 0.366762], - [0.944006, 0.377643, 0.365136], - [0.947180, 0.384178, 0.363701], - [0.950210, 0.390820, 0.362468], - [0.953099, 0.397563, 0.361438], - [0.955849, 0.404400, 0.360619], - [0.958464, 0.411324, 0.360014], - [0.960949, 0.418323, 0.359630], - [0.963310, 0.425390, 0.359469], - [0.965549, 0.432519, 0.359529], - [0.967671, 0.439703, 0.359810], - [0.969680, 0.446936, 0.360311], - [0.971582, 0.454210, 0.361030], - [0.973381, 0.461520, 0.361965], - [0.975082, 0.468861, 0.363111], - [0.976690, 0.476226, 0.364466], - [0.978210, 0.483612, 0.366025], - [0.979645, 0.491014, 0.367783], - [0.981000, 0.498428, 0.369734], - [0.982279, 0.505851, 0.371874], - [0.983485, 0.513280, 0.374198], - [0.984622, 0.520713, 0.376698], - [0.985693, 0.528148, 0.379371], - [0.986700, 0.535582, 0.382210], - [0.987646, 0.543015, 0.385210], - [0.988533, 0.550446, 0.388365], - [0.989363, 0.557873, 0.391671], - [0.990138, 0.565296, 0.395122], - [0.990871, 0.572706, 0.398714], - [0.991558, 0.580107, 0.402441], - [0.992196, 0.587502, 0.406299], - [0.992785, 0.594891, 0.410283], - [0.993326, 0.602275, 0.414390], - [0.993834, 0.609644, 0.418613], - [0.994309, 0.616999, 0.422950], - [0.994738, 0.624350, 0.427397], - [0.995122, 0.631696, 0.431951], - [0.995480, 0.639027, 0.436607], - [0.995810, 0.646344, 0.441361], - [0.996096, 0.653659, 0.446213], - [0.996341, 0.660969, 0.451160], - [0.996580, 0.668256, 0.456192], - [0.996775, 0.675541, 0.461314], - [0.996925, 0.682828, 0.466526], - [0.997077, 0.690088, 0.471811], - [0.997186, 0.697349, 0.477182], - [0.997254, 0.704611, 0.482635], - [0.997325, 0.711848, 0.488154], - [0.997351, 0.719089, 0.493755], - [0.997351, 0.726324, 0.499428], - [0.997341, 0.733545, 0.505167], - [0.997285, 0.740772, 0.510983], - [0.997228, 0.747981, 0.516859], - [0.997138, 0.755190, 0.522806], - [0.997019, 0.762398, 0.528821], - [0.996898, 0.769591, 0.534892], - [0.996727, 0.776795, 0.541039], - [0.996571, 0.783977, 0.547233], - [0.996369, 0.791167, 0.553499], - [0.996162, 0.798348, 0.559820], - [0.995932, 0.805527, 0.566202], - [0.995680, 0.812706, 0.572645], - [0.995424, 0.819875, 0.579140], - [0.995131, 0.827052, 0.585701], - [0.994851, 0.834213, 0.592307], - [0.994524, 0.841387, 0.598983], - [0.994222, 0.848540, 0.605696], - [0.993866, 0.855711, 0.612482], - [0.993545, 0.862859, 0.619299], - [0.993170, 0.870024, 0.626189], - [0.992831, 0.877168, 0.633109], - [0.992440, 0.884330, 0.640099], - [0.992089, 0.891470, 0.647116], - [0.991688, 0.898627, 0.654202], - [0.991332, 0.905763, 0.661309], - [0.990930, 0.912915, 0.668481], - [0.990570, 0.920049, 0.675675], - [0.990175, 0.927196, 0.682926], - [0.989815, 0.934329, 0.690198], - [0.989434, 0.941470, 0.697519], - [0.989077, 0.948604, 0.704863], - [0.988717, 0.955742, 0.712242], - [0.988367, 0.962878, 0.719649], - [0.988033, 0.970012, 0.727077], - [0.987691, 0.977154, 0.734536], - [0.987387, 0.984288, 0.742002], - [0.987053, 0.991438, 0.749504]] - -_inferno_data = [[0.001462, 0.000466, 0.013866], - [0.002267, 0.001270, 0.018570], - [0.003299, 0.002249, 0.024239], - [0.004547, 0.003392, 0.030909], - [0.006006, 0.004692, 0.038558], - [0.007676, 0.006136, 0.046836], - [0.009561, 0.007713, 0.055143], - [0.011663, 0.009417, 0.063460], - [0.013995, 0.011225, 0.071862], - [0.016561, 0.013136, 0.080282], - [0.019373, 0.015133, 0.088767], - [0.022447, 0.017199, 0.097327], - [0.025793, 0.019331, 0.105930], - [0.029432, 0.021503, 0.114621], - [0.033385, 0.023702, 0.123397], - [0.037668, 0.025921, 0.132232], - [0.042253, 0.028139, 0.141141], - [0.046915, 0.030324, 0.150164], - [0.051644, 0.032474, 0.159254], - [0.056449, 0.034569, 0.168414], - [0.061340, 0.036590, 0.177642], - [0.066331, 0.038504, 0.186962], - [0.071429, 0.040294, 0.196354], - [0.076637, 0.041905, 0.205799], - [0.081962, 0.043328, 0.215289], - [0.087411, 0.044556, 0.224813], - [0.092990, 0.045583, 0.234358], - [0.098702, 0.046402, 0.243904], - [0.104551, 0.047008, 0.253430], - [0.110536, 0.047399, 0.262912], - [0.116656, 0.047574, 0.272321], - [0.122908, 0.047536, 0.281624], - [0.129285, 0.047293, 0.290788], - [0.135778, 0.046856, 0.299776], - [0.142378, 0.046242, 0.308553], - [0.149073, 0.045468, 0.317085], - [0.155850, 0.044559, 0.325338], - [0.162689, 0.043554, 0.333277], - [0.169575, 0.042489, 0.340874], - [0.176493, 0.041402, 0.348111], - [0.183429, 0.040329, 0.354971], - [0.190367, 0.039309, 0.361447], - [0.197297, 0.038400, 0.367535], - [0.204209, 0.037632, 0.373238], - [0.211095, 0.037030, 0.378563], - [0.217949, 0.036615, 0.383522], - [0.224763, 0.036405, 0.388129], - [0.231538, 0.036405, 0.392400], - [0.238273, 0.036621, 0.396353], - [0.244967, 0.037055, 0.400007], - [0.251620, 0.037705, 0.403378], - [0.258234, 0.038571, 0.406485], - [0.264810, 0.039647, 0.409345], - [0.271347, 0.040922, 0.411976], - [0.277850, 0.042353, 0.414392], - [0.284321, 0.043933, 0.416608], - [0.290763, 0.045644, 0.418637], - [0.297178, 0.047470, 0.420491], - [0.303568, 0.049396, 0.422182], - [0.309935, 0.051407, 0.423721], - [0.316282, 0.053490, 0.425116], - [0.322610, 0.055634, 0.426377], - [0.328921, 0.057827, 0.427511], - [0.335217, 0.060060, 0.428524], - [0.341500, 0.062325, 0.429425], - [0.347771, 0.064616, 0.430217], - [0.354032, 0.066925, 0.430906], - [0.360284, 0.069247, 0.431497], - [0.366529, 0.071579, 0.431994], - [0.372768, 0.073915, 0.432400], - [0.379001, 0.076253, 0.432719], - [0.385228, 0.078591, 0.432955], - [0.391453, 0.080927, 0.433109], - [0.397674, 0.083257, 0.433183], - [0.403894, 0.085580, 0.433179], - [0.410113, 0.087896, 0.433098], - [0.416331, 0.090203, 0.432943], - [0.422549, 0.092501, 0.432714], - [0.428768, 0.094790, 0.432412], - [0.434987, 0.097069, 0.432039], - [0.441207, 0.099338, 0.431594], - [0.447428, 0.101597, 0.431080], - [0.453651, 0.103848, 0.430498], - [0.459875, 0.106089, 0.429846], - [0.466100, 0.108322, 0.429125], - [0.472328, 0.110547, 0.428334], - [0.478558, 0.112764, 0.427475], - [0.484789, 0.114974, 0.426548], - [0.491022, 0.117179, 0.425552], - [0.497257, 0.119379, 0.424488], - [0.503493, 0.121575, 0.423356], - [0.509730, 0.123769, 0.422156], - [0.515967, 0.125960, 0.420887], - [0.522206, 0.128150, 0.419549], - [0.528444, 0.130341, 0.418142], - [0.534683, 0.132534, 0.416667], - [0.540920, 0.134729, 0.415123], - [0.547157, 0.136929, 0.413511], - [0.553392, 0.139134, 0.411829], - [0.559624, 0.141346, 0.410078], - [0.565854, 0.143567, 0.408258], - [0.572081, 0.145797, 0.406369], - [0.578304, 0.148039, 0.404411], - [0.584521, 0.150294, 0.402385], - [0.590734, 0.152563, 0.400290], - [0.596940, 0.154848, 0.398125], - [0.603139, 0.157151, 0.395891], - [0.609330, 0.159474, 0.393589], - [0.615513, 0.161817, 0.391219], - [0.621685, 0.164184, 0.388781], - [0.627847, 0.166575, 0.386276], - [0.633998, 0.168992, 0.383704], - [0.640135, 0.171438, 0.381065], - [0.646260, 0.173914, 0.378359], - [0.652369, 0.176421, 0.375586], - [0.658463, 0.178962, 0.372748], - [0.664540, 0.181539, 0.369846], - [0.670599, 0.184153, 0.366879], - [0.676638, 0.186807, 0.363849], - [0.682656, 0.189501, 0.360757], - [0.688653, 0.192239, 0.357603], - [0.694627, 0.195021, 0.354388], - [0.700576, 0.197851, 0.351113], - [0.706500, 0.200728, 0.347777], - [0.712396, 0.203656, 0.344383], - [0.718264, 0.206636, 0.340931], - [0.724103, 0.209670, 0.337424], - [0.729909, 0.212759, 0.333861], - [0.735683, 0.215906, 0.330245], - [0.741423, 0.219112, 0.326576], - [0.747127, 0.222378, 0.322856], - [0.752794, 0.225706, 0.319085], - [0.758422, 0.229097, 0.315266], - [0.764010, 0.232554, 0.311399], - [0.769556, 0.236077, 0.307485], - [0.775059, 0.239667, 0.303526], - [0.780517, 0.243327, 0.299523], - [0.785929, 0.247056, 0.295477], - [0.791293, 0.250856, 0.291390], - [0.796607, 0.254728, 0.287264], - [0.801871, 0.258674, 0.283099], - [0.807082, 0.262692, 0.278898], - [0.812239, 0.266786, 0.274661], - [0.817341, 0.270954, 0.270390], - [0.822386, 0.275197, 0.266085], - [0.827372, 0.279517, 0.261750], - [0.832299, 0.283913, 0.257383], - [0.837165, 0.288385, 0.252988], - [0.841969, 0.292933, 0.248564], - [0.846709, 0.297559, 0.244113], - [0.851384, 0.302260, 0.239636], - [0.855992, 0.307038, 0.235133], - [0.860533, 0.311892, 0.230606], - [0.865006, 0.316822, 0.226055], - [0.869409, 0.321827, 0.221482], - [0.873741, 0.326906, 0.216886], - [0.878001, 0.332060, 0.212268], - [0.882188, 0.337287, 0.207628], - [0.886302, 0.342586, 0.202968], - [0.890341, 0.347957, 0.198286], - [0.894305, 0.353399, 0.193584], - [0.898192, 0.358911, 0.188860], - [0.902003, 0.364492, 0.184116], - [0.905735, 0.370140, 0.179350], - [0.909390, 0.375856, 0.174563], - [0.912966, 0.381636, 0.169755], - [0.916462, 0.387481, 0.164924], - [0.919879, 0.393389, 0.160070], - [0.923215, 0.399359, 0.155193], - [0.926470, 0.405389, 0.150292], - [0.929644, 0.411479, 0.145367], - [0.932737, 0.417627, 0.140417], - [0.935747, 0.423831, 0.135440], - [0.938675, 0.430091, 0.130438], - [0.941521, 0.436405, 0.125409], - [0.944285, 0.442772, 0.120354], - [0.946965, 0.449191, 0.115272], - [0.949562, 0.455660, 0.110164], - [0.952075, 0.462178, 0.105031], - [0.954506, 0.468744, 0.099874], - [0.956852, 0.475356, 0.094695], - [0.959114, 0.482014, 0.089499], - [0.961293, 0.488716, 0.084289], - [0.963387, 0.495462, 0.079073], - [0.965397, 0.502249, 0.073859], - [0.967322, 0.509078, 0.068659], - [0.969163, 0.515946, 0.063488], - [0.970919, 0.522853, 0.058367], - [0.972590, 0.529798, 0.053324], - [0.974176, 0.536780, 0.048392], - [0.975677, 0.543798, 0.043618], - [0.977092, 0.550850, 0.039050], - [0.978422, 0.557937, 0.034931], - [0.979666, 0.565057, 0.031409], - [0.980824, 0.572209, 0.028508], - [0.981895, 0.579392, 0.026250], - [0.982881, 0.586606, 0.024661], - [0.983779, 0.593849, 0.023770], - [0.984591, 0.601122, 0.023606], - [0.985315, 0.608422, 0.024202], - [0.985952, 0.615750, 0.025592], - [0.986502, 0.623105, 0.027814], - [0.986964, 0.630485, 0.030908], - [0.987337, 0.637890, 0.034916], - [0.987622, 0.645320, 0.039886], - [0.987819, 0.652773, 0.045581], - [0.987926, 0.660250, 0.051750], - [0.987945, 0.667748, 0.058329], - [0.987874, 0.675267, 0.065257], - [0.987714, 0.682807, 0.072489], - [0.987464, 0.690366, 0.079990], - [0.987124, 0.697944, 0.087731], - [0.986694, 0.705540, 0.095694], - [0.986175, 0.713153, 0.103863], - [0.985566, 0.720782, 0.112229], - [0.984865, 0.728427, 0.120785], - [0.984075, 0.736087, 0.129527], - [0.983196, 0.743758, 0.138453], - [0.982228, 0.751442, 0.147565], - [0.981173, 0.759135, 0.156863], - [0.980032, 0.766837, 0.166353], - [0.978806, 0.774545, 0.176037], - [0.977497, 0.782258, 0.185923], - [0.976108, 0.789974, 0.196018], - [0.974638, 0.797692, 0.206332], - [0.973088, 0.805409, 0.216877], - [0.971468, 0.813122, 0.227658], - [0.969783, 0.820825, 0.238686], - [0.968041, 0.828515, 0.249972], - [0.966243, 0.836191, 0.261534], - [0.964394, 0.843848, 0.273391], - [0.962517, 0.851476, 0.285546], - [0.960626, 0.859069, 0.298010], - [0.958720, 0.866624, 0.310820], - [0.956834, 0.874129, 0.323974], - [0.954997, 0.881569, 0.337475], - [0.953215, 0.888942, 0.351369], - [0.951546, 0.896226, 0.365627], - [0.950018, 0.903409, 0.380271], - [0.948683, 0.910473, 0.395289], - [0.947594, 0.917399, 0.410665], - [0.946809, 0.924168, 0.426373], - [0.946392, 0.930761, 0.442367], - [0.946403, 0.937159, 0.458592], - [0.946903, 0.943348, 0.474970], - [0.947937, 0.949318, 0.491426], - [0.949545, 0.955063, 0.507860], - [0.951740, 0.960587, 0.524203], - [0.954529, 0.965896, 0.540361], - [0.957896, 0.971003, 0.556275], - [0.961812, 0.975924, 0.571925], - [0.966249, 0.980678, 0.587206], - [0.971162, 0.985282, 0.602154], - [0.976511, 0.989753, 0.616760], - [0.982257, 0.994109, 0.631017], - [0.988362, 0.998364, 0.644924]] - -_plasma_data = [[0.050383, 0.029803, 0.527975], - [0.063536, 0.028426, 0.533124], - [0.075353, 0.027206, 0.538007], - [0.086222, 0.026125, 0.542658], - [0.096379, 0.025165, 0.547103], - [0.105980, 0.024309, 0.551368], - [0.115124, 0.023556, 0.555468], - [0.123903, 0.022878, 0.559423], - [0.132381, 0.022258, 0.563250], - [0.140603, 0.021687, 0.566959], - [0.148607, 0.021154, 0.570562], - [0.156421, 0.020651, 0.574065], - [0.164070, 0.020171, 0.577478], - [0.171574, 0.019706, 0.580806], - [0.178950, 0.019252, 0.584054], - [0.186213, 0.018803, 0.587228], - [0.193374, 0.018354, 0.590330], - [0.200445, 0.017902, 0.593364], - [0.207435, 0.017442, 0.596333], - [0.214350, 0.016973, 0.599239], - [0.221197, 0.016497, 0.602083], - [0.227983, 0.016007, 0.604867], - [0.234715, 0.015502, 0.607592], - [0.241396, 0.014979, 0.610259], - [0.248032, 0.014439, 0.612868], - [0.254627, 0.013882, 0.615419], - [0.261183, 0.013308, 0.617911], - [0.267703, 0.012716, 0.620346], - [0.274191, 0.012109, 0.622722], - [0.280648, 0.011488, 0.625038], - [0.287076, 0.010855, 0.627295], - [0.293478, 0.010213, 0.629490], - [0.299855, 0.009561, 0.631624], - [0.306210, 0.008902, 0.633694], - [0.312543, 0.008239, 0.635700], - [0.318856, 0.007576, 0.637640], - [0.325150, 0.006915, 0.639512], - [0.331426, 0.006261, 0.641316], - [0.337683, 0.005618, 0.643049], - [0.343925, 0.004991, 0.644710], - [0.350150, 0.004382, 0.646298], - [0.356359, 0.003798, 0.647810], - [0.362553, 0.003243, 0.649245], - [0.368733, 0.002724, 0.650601], - [0.374897, 0.002245, 0.651876], - [0.381047, 0.001814, 0.653068], - [0.387183, 0.001434, 0.654177], - [0.393304, 0.001114, 0.655199], - [0.399411, 0.000859, 0.656133], - [0.405503, 0.000678, 0.656977], - [0.411580, 0.000577, 0.657730], - [0.417642, 0.000564, 0.658390], - [0.423689, 0.000646, 0.658956], - [0.429719, 0.000831, 0.659425], - [0.435734, 0.001127, 0.659797], - [0.441732, 0.001540, 0.660069], - [0.447714, 0.002080, 0.660240], - [0.453677, 0.002755, 0.660310], - [0.459623, 0.003574, 0.660277], - [0.465550, 0.004545, 0.660139], - [0.471457, 0.005678, 0.659897], - [0.477344, 0.006980, 0.659549], - [0.483210, 0.008460, 0.659095], - [0.489055, 0.010127, 0.658534], - [0.494877, 0.011990, 0.657865], - [0.500678, 0.014055, 0.657088], - [0.506454, 0.016333, 0.656202], - [0.512206, 0.018833, 0.655209], - [0.517933, 0.021563, 0.654109], - [0.523633, 0.024532, 0.652901], - [0.529306, 0.027747, 0.651586], - [0.534952, 0.031217, 0.650165], - [0.540570, 0.034950, 0.648640], - [0.546157, 0.038954, 0.647010], - [0.551715, 0.043136, 0.645277], - [0.557243, 0.047331, 0.643443], - [0.562738, 0.051545, 0.641509], - [0.568201, 0.055778, 0.639477], - [0.573632, 0.060028, 0.637349], - [0.579029, 0.064296, 0.635126], - [0.584391, 0.068579, 0.632812], - [0.589719, 0.072878, 0.630408], - [0.595011, 0.077190, 0.627917], - [0.600266, 0.081516, 0.625342], - [0.605485, 0.085854, 0.622686], - [0.610667, 0.090204, 0.619951], - [0.615812, 0.094564, 0.617140], - [0.620919, 0.098934, 0.614257], - [0.625987, 0.103312, 0.611305], - [0.631017, 0.107699, 0.608287], - [0.636008, 0.112092, 0.605205], - [0.640959, 0.116492, 0.602065], - [0.645872, 0.120898, 0.598867], - [0.650746, 0.125309, 0.595617], - [0.655580, 0.129725, 0.592317], - [0.660374, 0.134144, 0.588971], - [0.665129, 0.138566, 0.585582], - [0.669845, 0.142992, 0.582154], - [0.674522, 0.147419, 0.578688], - [0.679160, 0.151848, 0.575189], - [0.683758, 0.156278, 0.571660], - [0.688318, 0.160709, 0.568103], - [0.692840, 0.165141, 0.564522], - [0.697324, 0.169573, 0.560919], - [0.701769, 0.174005, 0.557296], - [0.706178, 0.178437, 0.553657], - [0.710549, 0.182868, 0.550004], - [0.714883, 0.187299, 0.546338], - [0.719181, 0.191729, 0.542663], - [0.723444, 0.196158, 0.538981], - [0.727670, 0.200586, 0.535293], - [0.731862, 0.205013, 0.531601], - [0.736019, 0.209439, 0.527908], - [0.740143, 0.213864, 0.524216], - [0.744232, 0.218288, 0.520524], - [0.748289, 0.222711, 0.516834], - [0.752312, 0.227133, 0.513149], - [0.756304, 0.231555, 0.509468], - [0.760264, 0.235976, 0.505794], - [0.764193, 0.240396, 0.502126], - [0.768090, 0.244817, 0.498465], - [0.771958, 0.249237, 0.494813], - [0.775796, 0.253658, 0.491171], - [0.779604, 0.258078, 0.487539], - [0.783383, 0.262500, 0.483918], - [0.787133, 0.266922, 0.480307], - [0.790855, 0.271345, 0.476706], - [0.794549, 0.275770, 0.473117], - [0.798216, 0.280197, 0.469538], - [0.801855, 0.284626, 0.465971], - [0.805467, 0.289057, 0.462415], - [0.809052, 0.293491, 0.458870], - [0.812612, 0.297928, 0.455338], - [0.816144, 0.302368, 0.451816], - [0.819651, 0.306812, 0.448306], - [0.823132, 0.311261, 0.444806], - [0.826588, 0.315714, 0.441316], - [0.830018, 0.320172, 0.437836], - [0.833422, 0.324635, 0.434366], - [0.836801, 0.329105, 0.430905], - [0.840155, 0.333580, 0.427455], - [0.843484, 0.338062, 0.424013], - [0.846788, 0.342551, 0.420579], - [0.850066, 0.347048, 0.417153], - [0.853319, 0.351553, 0.413734], - [0.856547, 0.356066, 0.410322], - [0.859750, 0.360588, 0.406917], - [0.862927, 0.365119, 0.403519], - [0.866078, 0.369660, 0.400126], - [0.869203, 0.374212, 0.396738], - [0.872303, 0.378774, 0.393355], - [0.875376, 0.383347, 0.389976], - [0.878423, 0.387932, 0.386600], - [0.881443, 0.392529, 0.383229], - [0.884436, 0.397139, 0.379860], - [0.887402, 0.401762, 0.376494], - [0.890340, 0.406398, 0.373130], - [0.893250, 0.411048, 0.369768], - [0.896131, 0.415712, 0.366407], - [0.898984, 0.420392, 0.363047], - [0.901807, 0.425087, 0.359688], - [0.904601, 0.429797, 0.356329], - [0.907365, 0.434524, 0.352970], - [0.910098, 0.439268, 0.349610], - [0.912800, 0.444029, 0.346251], - [0.915471, 0.448807, 0.342890], - [0.918109, 0.453603, 0.339529], - [0.920714, 0.458417, 0.336166], - [0.923287, 0.463251, 0.332801], - [0.925825, 0.468103, 0.329435], - [0.928329, 0.472975, 0.326067], - [0.930798, 0.477867, 0.322697], - [0.933232, 0.482780, 0.319325], - [0.935630, 0.487712, 0.315952], - [0.937990, 0.492667, 0.312575], - [0.940313, 0.497642, 0.309197], - [0.942598, 0.502639, 0.305816], - [0.944844, 0.507658, 0.302433], - [0.947051, 0.512699, 0.299049], - [0.949217, 0.517763, 0.295662], - [0.951344, 0.522850, 0.292275], - [0.953428, 0.527960, 0.288883], - [0.955470, 0.533093, 0.285490], - [0.957469, 0.538250, 0.282096], - [0.959424, 0.543431, 0.278701], - [0.961336, 0.548636, 0.275305], - [0.963203, 0.553865, 0.271909], - [0.965024, 0.559118, 0.268513], - [0.966798, 0.564396, 0.265118], - [0.968526, 0.569700, 0.261721], - [0.970205, 0.575028, 0.258325], - [0.971835, 0.580382, 0.254931], - [0.973416, 0.585761, 0.251540], - [0.974947, 0.591165, 0.248151], - [0.976428, 0.596595, 0.244767], - [0.977856, 0.602051, 0.241387], - [0.979233, 0.607532, 0.238013], - [0.980556, 0.613039, 0.234646], - [0.981826, 0.618572, 0.231287], - [0.983041, 0.624131, 0.227937], - [0.984199, 0.629718, 0.224595], - [0.985301, 0.635330, 0.221265], - [0.986345, 0.640969, 0.217948], - [0.987332, 0.646633, 0.214648], - [0.988260, 0.652325, 0.211364], - [0.989128, 0.658043, 0.208100], - [0.989935, 0.663787, 0.204859], - [0.990681, 0.669558, 0.201642], - [0.991365, 0.675355, 0.198453], - [0.991985, 0.681179, 0.195295], - [0.992541, 0.687030, 0.192170], - [0.993032, 0.692907, 0.189084], - [0.993456, 0.698810, 0.186041], - [0.993814, 0.704741, 0.183043], - [0.994103, 0.710698, 0.180097], - [0.994324, 0.716681, 0.177208], - [0.994474, 0.722691, 0.174381], - [0.994553, 0.728728, 0.171622], - [0.994561, 0.734791, 0.168938], - [0.994495, 0.740880, 0.166335], - [0.994355, 0.746995, 0.163821], - [0.994141, 0.753137, 0.161404], - [0.993851, 0.759304, 0.159092], - [0.993482, 0.765499, 0.156891], - [0.993033, 0.771720, 0.154808], - [0.992505, 0.777967, 0.152855], - [0.991897, 0.784239, 0.151042], - [0.991209, 0.790537, 0.149377], - [0.990439, 0.796859, 0.147870], - [0.989587, 0.803205, 0.146529], - [0.988648, 0.809579, 0.145357], - [0.987621, 0.815978, 0.144363], - [0.986509, 0.822401, 0.143557], - [0.985314, 0.828846, 0.142945], - [0.984031, 0.835315, 0.142528], - [0.982653, 0.841812, 0.142303], - [0.981190, 0.848329, 0.142279], - [0.979644, 0.854866, 0.142453], - [0.977995, 0.861432, 0.142808], - [0.976265, 0.868016, 0.143351], - [0.974443, 0.874622, 0.144061], - [0.972530, 0.881250, 0.144923], - [0.970533, 0.887896, 0.145919], - [0.968443, 0.894564, 0.147014], - [0.966271, 0.901249, 0.148180], - [0.964021, 0.907950, 0.149370], - [0.961681, 0.914672, 0.150520], - [0.959276, 0.921407, 0.151566], - [0.956808, 0.928152, 0.152409], - [0.954287, 0.934908, 0.152921], - [0.951726, 0.941671, 0.152925], - [0.949151, 0.948435, 0.152178], - [0.946602, 0.955190, 0.150328], - [0.944152, 0.961916, 0.146861], - [0.941896, 0.968590, 0.140956], - [0.940015, 0.975158, 0.131326]] - -_viridis_data = [[0.267004, 0.004874, 0.329415], - [0.268510, 0.009605, 0.335427], - [0.269944, 0.014625, 0.341379], - [0.271305, 0.019942, 0.347269], - [0.272594, 0.025563, 0.353093], - [0.273809, 0.031497, 0.358853], - [0.274952, 0.037752, 0.364543], - [0.276022, 0.044167, 0.370164], - [0.277018, 0.050344, 0.375715], - [0.277941, 0.056324, 0.381191], - [0.278791, 0.062145, 0.386592], - [0.279566, 0.067836, 0.391917], - [0.280267, 0.073417, 0.397163], - [0.280894, 0.078907, 0.402329], - [0.281446, 0.084320, 0.407414], - [0.281924, 0.089666, 0.412415], - [0.282327, 0.094955, 0.417331], - [0.282656, 0.100196, 0.422160], - [0.282910, 0.105393, 0.426902], - [0.283091, 0.110553, 0.431554], - [0.283197, 0.115680, 0.436115], - [0.283229, 0.120777, 0.440584], - [0.283187, 0.125848, 0.444960], - [0.283072, 0.130895, 0.449241], - [0.282884, 0.135920, 0.453427], - [0.282623, 0.140926, 0.457517], - [0.282290, 0.145912, 0.461510], - [0.281887, 0.150881, 0.465405], - [0.281412, 0.155834, 0.469201], - [0.280868, 0.160771, 0.472899], - [0.280255, 0.165693, 0.476498], - [0.279574, 0.170599, 0.479997], - [0.278826, 0.175490, 0.483397], - [0.278012, 0.180367, 0.486697], - [0.277134, 0.185228, 0.489898], - [0.276194, 0.190074, 0.493001], - [0.275191, 0.194905, 0.496005], - [0.274128, 0.199721, 0.498911], - [0.273006, 0.204520, 0.501721], - [0.271828, 0.209303, 0.504434], - [0.270595, 0.214069, 0.507052], - [0.269308, 0.218818, 0.509577], - [0.267968, 0.223549, 0.512008], - [0.266580, 0.228262, 0.514349], - [0.265145, 0.232956, 0.516599], - [0.263663, 0.237631, 0.518762], - [0.262138, 0.242286, 0.520837], - [0.260571, 0.246922, 0.522828], - [0.258965, 0.251537, 0.524736], - [0.257322, 0.256130, 0.526563], - [0.255645, 0.260703, 0.528312], - [0.253935, 0.265254, 0.529983], - [0.252194, 0.269783, 0.531579], - [0.250425, 0.274290, 0.533103], - [0.248629, 0.278775, 0.534556], - [0.246811, 0.283237, 0.535941], - [0.244972, 0.287675, 0.537260], - [0.243113, 0.292092, 0.538516], - [0.241237, 0.296485, 0.539709], - [0.239346, 0.300855, 0.540844], - [0.237441, 0.305202, 0.541921], - [0.235526, 0.309527, 0.542944], - [0.233603, 0.313828, 0.543914], - [0.231674, 0.318106, 0.544834], - [0.229739, 0.322361, 0.545706], - [0.227802, 0.326594, 0.546532], - [0.225863, 0.330805, 0.547314], - [0.223925, 0.334994, 0.548053], - [0.221989, 0.339161, 0.548752], - [0.220057, 0.343307, 0.549413], - [0.218130, 0.347432, 0.550038], - [0.216210, 0.351535, 0.550627], - [0.214298, 0.355619, 0.551184], - [0.212395, 0.359683, 0.551710], - [0.210503, 0.363727, 0.552206], - [0.208623, 0.367752, 0.552675], - [0.206756, 0.371758, 0.553117], - [0.204903, 0.375746, 0.553533], - [0.203063, 0.379716, 0.553925], - [0.201239, 0.383670, 0.554294], - [0.199430, 0.387607, 0.554642], - [0.197636, 0.391528, 0.554969], - [0.195860, 0.395433, 0.555276], - [0.194100, 0.399323, 0.555565], - [0.192357, 0.403199, 0.555836], - [0.190631, 0.407061, 0.556089], - [0.188923, 0.410910, 0.556326], - [0.187231, 0.414746, 0.556547], - [0.185556, 0.418570, 0.556753], - [0.183898, 0.422383, 0.556944], - [0.182256, 0.426184, 0.557120], - [0.180629, 0.429975, 0.557282], - [0.179019, 0.433756, 0.557430], - [0.177423, 0.437527, 0.557565], - [0.175841, 0.441290, 0.557685], - [0.174274, 0.445044, 0.557792], - [0.172719, 0.448791, 0.557885], - [0.171176, 0.452530, 0.557965], - [0.169646, 0.456262, 0.558030], - [0.168126, 0.459988, 0.558082], - [0.166617, 0.463708, 0.558119], - [0.165117, 0.467423, 0.558141], - [0.163625, 0.471133, 0.558148], - [0.162142, 0.474838, 0.558140], - [0.160665, 0.478540, 0.558115], - [0.159194, 0.482237, 0.558073], - [0.157729, 0.485932, 0.558013], - [0.156270, 0.489624, 0.557936], - [0.154815, 0.493313, 0.557840], - [0.153364, 0.497000, 0.557724], - [0.151918, 0.500685, 0.557587], - [0.150476, 0.504369, 0.557430], - [0.149039, 0.508051, 0.557250], - [0.147607, 0.511733, 0.557049], - [0.146180, 0.515413, 0.556823], - [0.144759, 0.519093, 0.556572], - [0.143343, 0.522773, 0.556295], - [0.141935, 0.526453, 0.555991], - [0.140536, 0.530132, 0.555659], - [0.139147, 0.533812, 0.555298], - [0.137770, 0.537492, 0.554906], - [0.136408, 0.541173, 0.554483], - [0.135066, 0.544853, 0.554029], - [0.133743, 0.548535, 0.553541], - [0.132444, 0.552216, 0.553018], - [0.131172, 0.555899, 0.552459], - [0.129933, 0.559582, 0.551864], - [0.128729, 0.563265, 0.551229], - [0.127568, 0.566949, 0.550556], - [0.126453, 0.570633, 0.549841], - [0.125394, 0.574318, 0.549086], - [0.124395, 0.578002, 0.548287], - [0.123463, 0.581687, 0.547445], - [0.122606, 0.585371, 0.546557], - [0.121831, 0.589055, 0.545623], - [0.121148, 0.592739, 0.544641], - [0.120565, 0.596422, 0.543611], - [0.120092, 0.600104, 0.542530], - [0.119738, 0.603785, 0.541400], - [0.119512, 0.607464, 0.540218], - [0.119423, 0.611141, 0.538982], - [0.119483, 0.614817, 0.537692], - [0.119699, 0.618490, 0.536347], - [0.120081, 0.622161, 0.534946], - [0.120638, 0.625828, 0.533488], - [0.121380, 0.629492, 0.531973], - [0.122312, 0.633153, 0.530398], - [0.123444, 0.636809, 0.528763], - [0.124780, 0.640461, 0.527068], - [0.126326, 0.644107, 0.525311], - [0.128087, 0.647749, 0.523491], - [0.130067, 0.651384, 0.521608], - [0.132268, 0.655014, 0.519661], - [0.134692, 0.658636, 0.517649], - [0.137339, 0.662252, 0.515571], - [0.140210, 0.665859, 0.513427], - [0.143303, 0.669459, 0.511215], - [0.146616, 0.673050, 0.508936], - [0.150148, 0.676631, 0.506589], - [0.153894, 0.680203, 0.504172], - [0.157851, 0.683765, 0.501686], - [0.162016, 0.687316, 0.499129], - [0.166383, 0.690856, 0.496502], - [0.170948, 0.694384, 0.493803], - [0.175707, 0.697900, 0.491033], - [0.180653, 0.701402, 0.488189], - [0.185783, 0.704891, 0.485273], - [0.191090, 0.708366, 0.482284], - [0.196571, 0.711827, 0.479221], - [0.202219, 0.715272, 0.476084], - [0.208030, 0.718701, 0.472873], - [0.214000, 0.722114, 0.469588], - [0.220124, 0.725509, 0.466226], - [0.226397, 0.728888, 0.462789], - [0.232815, 0.732247, 0.459277], - [0.239374, 0.735588, 0.455688], - [0.246070, 0.738910, 0.452024], - [0.252899, 0.742211, 0.448284], - [0.259857, 0.745492, 0.444467], - [0.266941, 0.748751, 0.440573], - [0.274149, 0.751988, 0.436601], - [0.281477, 0.755203, 0.432552], - [0.288921, 0.758394, 0.428426], - [0.296479, 0.761561, 0.424223], - [0.304148, 0.764704, 0.419943], - [0.311925, 0.767822, 0.415586], - [0.319809, 0.770914, 0.411152], - [0.327796, 0.773980, 0.406640], - [0.335885, 0.777018, 0.402049], - [0.344074, 0.780029, 0.397381], - [0.352360, 0.783011, 0.392636], - [0.360741, 0.785964, 0.387814], - [0.369214, 0.788888, 0.382914], - [0.377779, 0.791781, 0.377939], - [0.386433, 0.794644, 0.372886], - [0.395174, 0.797475, 0.367757], - [0.404001, 0.800275, 0.362552], - [0.412913, 0.803041, 0.357269], - [0.421908, 0.805774, 0.351910], - [0.430983, 0.808473, 0.346476], - [0.440137, 0.811138, 0.340967], - [0.449368, 0.813768, 0.335384], - [0.458674, 0.816363, 0.329727], - [0.468053, 0.818921, 0.323998], - [0.477504, 0.821444, 0.318195], - [0.487026, 0.823929, 0.312321], - [0.496615, 0.826376, 0.306377], - [0.506271, 0.828786, 0.300362], - [0.515992, 0.831158, 0.294279], - [0.525776, 0.833491, 0.288127], - [0.535621, 0.835785, 0.281908], - [0.545524, 0.838039, 0.275626], - [0.555484, 0.840254, 0.269281], - [0.565498, 0.842430, 0.262877], - [0.575563, 0.844566, 0.256415], - [0.585678, 0.846661, 0.249897], - [0.595839, 0.848717, 0.243329], - [0.606045, 0.850733, 0.236712], - [0.616293, 0.852709, 0.230052], - [0.626579, 0.854645, 0.223353], - [0.636902, 0.856542, 0.216620], - [0.647257, 0.858400, 0.209861], - [0.657642, 0.860219, 0.203082], - [0.668054, 0.861999, 0.196293], - [0.678489, 0.863742, 0.189503], - [0.688944, 0.865448, 0.182725], - [0.699415, 0.867117, 0.175971], - [0.709898, 0.868751, 0.169257], - [0.720391, 0.870350, 0.162603], - [0.730889, 0.871916, 0.156029], - [0.741388, 0.873449, 0.149561], - [0.751884, 0.874951, 0.143228], - [0.762373, 0.876424, 0.137064], - [0.772852, 0.877868, 0.131109], - [0.783315, 0.879285, 0.125405], - [0.793760, 0.880678, 0.120005], - [0.804182, 0.882046, 0.114965], - [0.814576, 0.883393, 0.110347], - [0.824940, 0.884720, 0.106217], - [0.835270, 0.886029, 0.102646], - [0.845561, 0.887322, 0.099702], - [0.855810, 0.888601, 0.097452], - [0.866013, 0.889868, 0.095953], - [0.876168, 0.891125, 0.095250], - [0.886271, 0.892374, 0.095374], - [0.896320, 0.893616, 0.096335], - [0.906311, 0.894855, 0.098125], - [0.916242, 0.896091, 0.100717], - [0.926106, 0.897330, 0.104071], - [0.935904, 0.898570, 0.108131], - [0.945636, 0.899815, 0.112838], - [0.955300, 0.901065, 0.118128], - [0.964894, 0.902323, 0.123941], - [0.974417, 0.903590, 0.130215], - [0.983868, 0.904867, 0.136897], - [0.993248, 0.906157, 0.143936]] - -from matplotlib.colors import ListedColormap - -cmaps = {} -for (name, data) in (('magma', _magma_data), - ('inferno', _inferno_data), - ('plasma', _plasma_data), - ('viridis', _viridis_data)): - - cmaps[name] = ListedColormap(data, name=name) - -magma = cmaps['magma'] -inferno = cmaps['inferno'] -plasma = cmaps['plasma'] -viridis = cmaps['viridis'] diff --git a/common_grid.py b/common_grid.py index be4b9b8..ead74af 100644 --- a/common_grid.py +++ b/common_grid.py @@ -6,6 +6,15 @@ from monthly_avg_roms import * from monthly_avg_cice import * +# Interpolate ROMS output to a regular quarter-degree grid for easy comparison +# with FESOM. Write monthly averages of surface heat and salt flux, the +# surface stress vector and its curl, and the sea ice velocity vector. +# Input: +# roms_file = path to ROMS output file containing 5-day averages for the entire +# simulation, starting on 1 January +# cice_file = path to CICE output file containing 5-day averages for the same +# timesteps as ROMS_file +# out_file = path to desired output file def common_grid (roms_file, cice_file, out_file): # Resolution of common grid (degrees, same for lat and lon) diff --git a/compare_nic_seasonal.py b/compare_nic_seasonal.py index 1707f74..815e06b 100644 --- a/compare_nic_seasonal.py +++ b/compare_nic_seasonal.py @@ -4,20 +4,40 @@ from rotate_vector_cice import * from seasonal_avg_cice import * +# Make a 4x2 plot comparing seasonal averages of the given variable from my +# implementation of CICE in MetROMS (bottom) and Nic Hannah's implementation +# coupled to MOM (top). +# Input: +# cice_file = path to CICE file from MetROMS simulation, containing at least +# one complete December-November period, in 5-day averages. If +# there are multiple such periods, the last one will be used. +# var_name = variable name to plot +# colour_bounds = optional array of size 2 containing bounds on colour scale +# save = optional boolean indicating to save the file rather than display it +# on screen +# fig_name = if save=True, filename for figure def compare_nic_seasonal (cice_file, var_name, colour_bounds=None, save=False, fig_name=None): + # Piece together the paths to Nic's monthly averaged output on raijin nic_dir_head = '/g/data/gh5/access_om_025-CORE_NYF/output' output_number = 137 nic_dir_tail = '/ice/HISTORY/' nic_file_head = 'iceh.0' nic_year_number = 133 + # Maximum j-index to read in Nic's output max_j = 300 + # Look at the variable in my output id = Dataset(cice_file, 'r') + # Save units units = id.variables[var_name].units + # Check if this is a vector we need to rotate to lon-lat space if var_name in ['uvel', 'vvel', 'strairx', 'strairy', 'strocnx', 'strocny']: rotate = True + # Read the angle of grid rotation angle = id.variables['ANGLE'][:-15,:] + # Figure out whether this is an x or y component, and what the name of + # the other component is if var_name == 'uvel': cmp_flag = 'x' other_name = 'vvel' @@ -32,6 +52,7 @@ def compare_nic_seasonal (cice_file, var_name, colour_bounds=None, save=False, f other_name = var_name.replace('x', 'y') else: rotate = False + # Read the correct grid (tracer or velocity) grid_string = id.variables[var_name].coordinates if grid_string.startswith('ULON'): lon_name = 'ULON' @@ -42,13 +63,14 @@ def compare_nic_seasonal (cice_file, var_name, colour_bounds=None, save=False, f id.close() # Number of days in each month (this is just for Nic's output) + # Note Nic doesn't run with leap years ndays_month = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31] # Season names for titles season_names = ['DJF', 'MAM', 'JJA', 'SON'] # Degrees to radians conversion deg2rad = pi/180.0 - # Read the CICE grid + # Read my CICE grid id = Dataset(cice_file, 'r') cice_lon_tmp = id.variables[lon_name][:-15,:] cice_lat_tmp = id.variables[lat_name][:-15,:] @@ -65,8 +87,10 @@ def compare_nic_seasonal (cice_file, var_name, colour_bounds=None, save=False, f # Get seasonal averages of CICE data if rotate: + # Average both components of the vector this_cmp = seasonal_avg_cice(cice_file, var_name, [num_lat, num_lon]) other_cmp = seasonal_avg_cice(cice_file, other_name, [num_lat, num_lon]) + # Rotate to lon-lat space if cmp_flag == 'x': cice_data_tmp, other_tmp = rotate_vector_cice(this_cmp, other_cmp, angle) elif cmp_flag == 'y': @@ -87,9 +111,12 @@ def compare_nic_seasonal (cice_file, var_name, colour_bounds=None, save=False, f nic_lat = id.variables[lat_name][:max_j,:] id.close() + # Get seasonal averages of Nic's output nic_data = ma.empty([4, size(nic_lon,0), size(nic_lon,1)]) nic_data[:,:,:] = 0.0 + # Loop over seasons for season in range(4): + # Figure out what months we care about for this season if season == 0: months = [12, 1, 2] elif season == 1: @@ -98,10 +125,12 @@ def compare_nic_seasonal (cice_file, var_name, colour_bounds=None, save=False, f months = [6, 7, 8] elif season == 3: months = [9, 10, 11] + # Days in season so far season_days = 0 - + # Loop over months for month in months: if month == 12: + # Read December from the previous year filename = nic_dir_head + str(output_number-1) + nic_dir_tail + nic_file_head + str(nic_year_number-1) + '-' + str(month) + '.nc' else: if month < 10: @@ -109,8 +138,10 @@ def compare_nic_seasonal (cice_file, var_name, colour_bounds=None, save=False, f else: filename = nic_dir_head + str(output_number) + nic_dir_tail + nic_file_head + str(nic_year_number) + '-' + str(month) + '.nc' id = Dataset(filename, 'r') + # Integrate over time nic_data[season,:,:] += id.variables[var_name][0,:max_j,:]*ndays_month[month-1] season_days += ndays_month[month-1] + # Convert from integral to average nic_data[season,:,:] /= season_days # Convert both grids to spherical coordinates @@ -119,30 +150,40 @@ def compare_nic_seasonal (cice_file, var_name, colour_bounds=None, save=False, f nic_x = -(nic_lat+90)*cos(nic_lon*deg2rad+pi/2) nic_y = (nic_lat+90)*sin(nic_lon*deg2rad+pi/2) + # Boundaries on plots bdry1 = -35 bdry2 = 39 bdry3 = -35 bdry4 = 39 if colour_bounds is not None: + # User-defined colour bounds lev = linspace(colour_bounds[0], colour_bounds[1], num=50) if colour_bounds[0] == -colour_bounds[1]: + # Centered on zero; use red-yellow-blue colourmap colour_map = 'RdYlBu_r' else: + # Not centered on zero; go for a rainbow colour_map = 'jet' else: + # Automatic colour bounds based on min/max of the data if var_name in ['uvel', 'vvel', 'strairx', 'strairy', 'strocnx', 'strocny']: + # Center on zero and use red-yellow-blue colourmap max_val = max(amax(abs(nic_data)), amax(abs(cice_data))) lev = linspace(-max_val, max_val, num=50) colour_map = 'RdYlBu_r' else: + # Not centered on zero; go for a rainbow min_val = min(amin(nic_data), amin(cice_data)) max_val = max(amax(nic_data), amax(cice_data)) lev = linspace(min_val, max_val, num=50) colour_map = 'jet' + # Make the figure fig = figure(figsize=(20,9)) + # Loop over seasons for season in range(4): + # Nic's output ax = fig.add_subplot(2, 4, season+1, aspect='equal') contourf(nic_x, nic_y, nic_data[season,:,:], lev, cmap=colour_map, extend='both') if season == 0: @@ -151,6 +192,7 @@ def compare_nic_seasonal (cice_file, var_name, colour_bounds=None, save=False, f xlim([bdry1, bdry2]) ylim([bdry3, bdry4]) axis('off') + # My output ax = fig.add_subplot(2, 4, season+5, aspect='equal') img = contourf(cice_x, cice_y, cice_data[season,:,:], lev, cmap=colour_map, extend='both') if season == 0: @@ -158,9 +200,11 @@ def compare_nic_seasonal (cice_file, var_name, colour_bounds=None, save=False, f xlim([bdry1, bdry2]) ylim([bdry3, bdry4]) axis('off') + # Add colourbar on the bottom cbaxes = fig.add_axes([0.25, 0.04, 0.5, 0.02]) cbar = colorbar(img, orientation='horizontal', cax=cbaxes) cbar.ax.tick_params(labelsize=16) + # Main title with variable name and units suptitle(var_name + ' (' + units + ')', fontsize=30) subplots_adjust(wspace=0.025,hspace=0.025) @@ -171,6 +215,7 @@ def compare_nic_seasonal (cice_file, var_name, colour_bounds=None, save=False, f fig.show() +# Command-line interface if __name__ == "__main__": cice_file = raw_input("Path to CICE file, containing at least one complete Dec-Nov period: ") @@ -188,14 +233,18 @@ def compare_nic_seasonal (cice_file, var_name, colour_bounds=None, save=False, f elif action == 'd': save = False fig_name = None + # Make the plot compare_nic_seasonal(cice_file, var_name, colour_bounds, save, fig_name) + # Repeat until the user is finished while True: repeat = raw_input("Make another plot (y/n)? ") if repeat == 'y': + # Ask for changes to input parameters until the user is finished while True: changes = raw_input("Enter a parameter to change: (1) file path, (2) variable name, (3) colour bounds, (4) save/display; or enter to continue: ") if len(changes) == 0: + # No more changes to parameters break else: if int(changes) == 1: @@ -212,9 +261,12 @@ def compare_nic_seasonal (cice_file, var_name, colour_bounds=None, save=False, f elif int(changes) == 4: save = not save if save: + # Get new figure name fig_name = raw_input("File name for figure: ") + # Make the plot compare_nic_seasonal(cice_file, var_name, colour_bounds, save, fig_name) else: + # No more plots break diff --git a/convert_era.job b/convert_era.job index fd9261d..3a557c2 100644 --- a/convert_era.job +++ b/convert_era.job @@ -7,7 +7,7 @@ #PBS -v YEAR,COUNT # Call romscice_atm_subdaily.py for the given year, and self-submit enough times -# to process the entire ECCO2 file for that year (the python script only +# to process the entire ERA-Interim file for that year (the python script only # processes 100 6-hour timesteps at once, otherwise memory overflows). # To get it started for eg 1992, type diff --git a/era_evap.job b/era_evap.job index dc42ca7..4306fd4 100644 --- a/era_evap.job +++ b/era_evap.job @@ -6,6 +6,10 @@ #PBS -j oe #PBS -v YEAR,COUNT +# Call romscice_evap.py for the given year, and self-submit enough times to +# process the entire ERA-Interim file for that year (the python script only +# process 50 12-hour timesteps at once, otherwise memory overflows). + # To get it started for eg 1992, type # qsub -v YEAR=1992,COUNT=0 era_evap.job diff --git a/era_gpcp.py b/era_gpcp.py index 37fbf4e..ede6a22 100644 --- a/era_gpcp.py +++ b/era_gpcp.py @@ -2,30 +2,40 @@ from netCDF4 import Dataset from matplotlib.pyplot import * +# Plot precipitation (1992-2005 average) in the ERA-Interim versus GPCP datasets +# over the Southern Ocean. def era_gpcp (): + # Beginning and end of paths to ERA-Interim precip files, monthly averaged, + # on original grid era_head = '/short/y99/kaa561/CMIP5_forcing/atmos/climatology/ERA_Interim_monthly/FC_' era_tail = '_monthly_orig.nc' + # Beginning of paths to GPCP files gpcp_head = '/short/m68/kaa561/gpcp/gpcp_cdr_v23rB1_y' + # Days per month (leap years will be handled later) days_per_month = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31] start_year = 1992 end_year = 2005 deg2rad = pi/180.0 + # Read each grid id = Dataset(gpcp_head + str(start_year) + '_m01.nc', 'r') gpcp_lat_1d = id.variables['latitude'][:] gpcp_lon_1d = id.variables['longitude'][:] id.close() - id = Dataset(era_head + str(start_year) + era_tail, 'r') era_lat_1d = id.variables['latitude'][:] era_lon_1d = id.variables['longitude'][:] id.close() + # Start integrating precipitation data gpcp_precip = zeros([size(gpcp_lat_1d), size(gpcp_lon_1d)]) era_precip = zeros([size(era_lat_1d), size(era_lon_1d)]) + # Number of days so far ndays = 0 + # Loop over years for year in range(start_year, end_year+1): + # Check if it's a leap year leap_year = False if mod(year, 4) == 0: leap_year = True @@ -37,54 +47,69 @@ def era_gpcp (): days_per_month[1] = 29 else: days_per_month[1] = 28 + # Loop over months for month in range(12): if month + 1 < 10: month_string = '0' + str(month+1) else: month_string = str(month+1) + # Read GPCP data id = Dataset(gpcp_head + str(year) + '_m' + month_string + '.nc', 'r') gpcp_precip += id.variables['precip'][:,:]*days_per_month[month] id.close() + # Read ERA-Interim data id = Dataset(era_head + str(year) + era_tail, 'r') era_precip += id.variables['tp'][month,:,:]*days_per_month[month] id.close() ndays += days_per_month[month] + # Convert from integral to average gpcp_precip /= ndays era_precip /= ndays - gpcp_precip *= 1e-3*365.25 - era_precip *= 2*365.25 + # Convert to m/y + gpcp_precip *= 1e-3*365.25 # originally mm/day + era_precip *= 2*365.25 # originally m/12h + # Make 2D versions of lat and lon arrays gpcp_lon, gpcp_lat = meshgrid(gpcp_lon_1d, gpcp_lat_1d) era_lon, era_lat = meshgrid(era_lon_1d, era_lat_1d) + # Polar coordinates for plotting gpcp_x = -(gpcp_lat+90)*cos(gpcp_lon*deg2rad+pi/2) gpcp_y = (gpcp_lat+90)*sin(gpcp_lon*deg2rad+pi/2) era_x = -(era_lat+90)*cos(era_lon*deg2rad+pi/2) era_y = (era_lat+90)*sin(era_lon*deg2rad+pi/2) + # Colour levels lev = linspace(0, 2, num=50) + # Northern boundary 40S bdry = -40+90 + # Make figure fig = figure(figsize=(20,9)) + # GPCP fig.add_subplot(1,2,1, aspect='equal') contourf(gpcp_x, gpcp_y, gpcp_precip, lev, extend='max') title('GPCP', fontsize=24) xlim([-bdry, bdry]) ylim([-bdry, bdry]) axis('off') + # ERA-Interim fig.add_subplot(1,2,2, aspect='equal') img = contourf(era_x, era_y, era_precip, lev, extend='max') title('ERA-Interim', fontsize=24) xlim([-bdry, bdry]) ylim([-bdry, bdry]) axis('off') + # Add colourbar on the bottom cbaxes = fig.add_axes([0.3, 0.04, 0.4, 0.04]) cbar = colorbar(img, orientation='horizontal', cax=cbaxes) + # Main title suptitle('Precipitation (m/y), 1992-2005 average', fontsize=30) fig.savefig('era_gpcp.png') +# Command-line interface if __name__ == "__main__": era_gpcp() diff --git a/evap_era_monthly.py b/evap_era_monthly.py index 93cb288..e22bb94 100644 --- a/evap_era_monthly.py +++ b/evap_era_monthly.py @@ -2,8 +2,25 @@ from netCDF4 import Dataset, num2date from matplotlib.pyplot import * -def evap_era_monthly (roms_file, month, colour_bounds=None, save=False, fig_name=None): +# Plot monthly-averaged evaporation as calculated by MetROMS (ROMS and CICE +# calculate separately, so merge them based on sea ice concentration) and +# compare it to ERA-Interim for the given month. +# Input: +# roms_file = path to ROMS output file containing 5-day averages. If more than +# one instance of the given month is present, the last one will +# be used. Note that the evaporation variable will only be output +# if the EMINUSP option is on. +# cice_file = path to CICE output file, containing 5-day averages for the same +# timesteps as ROMS. +# month = month to plot (0 to 11) +# colour_bounds = optional array of size 2 containing bounds on colour scale +# save = optional boolean indicating to save the figure rather than display it +# on screen +# fig_name = if save=True, filename for figure +def evap_era_monthly (roms_file, cice_file, month, colour_bounds=None, save=False, fig_name=None): + # Beginning and end of paths to ERA-Interim monthly averaged files + # containing evaporation era_head = '/short/y99/kaa561/CMIP5_forcing/atmos/climatology/ERA_Interim_monthly/ER_' era_tail = '_monthly_orig.nc' # Starting and ending days in each month @@ -14,21 +31,29 @@ def evap_era_monthly (roms_file, month, colour_bounds=None, save=False, fig_name month_name = ['January', 'February', 'March', 'April', 'May', 'June', 'July', 'August', 'September', 'October', 'November', 'December'] # Degrees to radians conversion deg2rad = pi/180.0 - # Conversion factor: m/12h to kg/m^2/2, opposite signs + # Conversion factor: m/12h to kg/m^2/s, opposite signs era_conv = -1000./(12.*60.*60.) + # Conversion factor: cm/day to kg/m^2/s, opposite signs + cice_conv = -1000/(100.*24.*60.*60.) - id = Dataset(roms_file, 'r') - roms_lon = id.variables['lon_rho'][:-15,:-1] - roms_lat = id.variables['lat_rho'][:-15,:-1] - time_id = id.variables['ocean_time'] + # Read ROMS grid and time axis + id_roms = Dataset(roms_file, 'r') + roms_lon = id_roms.variables['lon_rho'][1:-16,1:-1] + roms_lat = id_roms.variables['lat_rho'][1:-16,1:-1] + time_id = id_roms.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 time = num2date(time_id[:], units=time_id.units, calendar=time_id.calendar.lower()) + # Open the CICE file + id_cice = Dataset(cice_file, 'r') + next_month = mod(month+1,12) prev_month = mod(month-1,12) - end_t = -1 + # Loop backwards from the end of the file to find the last timestep we + # care about (end of the month for the last instance of this month) + end_t = -1 # Missing value flag for t in range(size(time)-1, -1, -1): if time[t].month-1 == month and time[t].day in range(end_day[month]-2, end_day[month]+1): end_t = t @@ -36,11 +61,14 @@ def evap_era_monthly (roms_file, month, colour_bounds=None, save=False, fig_name if time[t].month-1 == next_month and time[t].day in range(start_day[next_month], start_day[next_month]+2): end_t = t break + # Check if we actually found it if end_t == -1: print 'Error: ' + roms_file + ' does not contain a complete ' + month_name[month] return - start_t = -1 + # Now loop backwards to find the first timestep we care about (beginning of + # the same month) + start_t = -1 # Missing value flag for t in range(end_t, -1, -1): if time[t].month-1 == prev_month and time[t].day in range(end_day[prev_month]-1, end_day[prev_month]+1): start_t = t @@ -48,14 +76,21 @@ def evap_era_monthly (roms_file, month, colour_bounds=None, save=False, fig_name if time[t].month-1 == month and time[t].day in range(start_day[month], start_day[month]+3): start_t = t break + # Check if we actually found it if start_t == -1: print 'Error: ' + roms_file + ' does not contain a complete ' + month_name[month] return - leap_year = False + # Figure out what year this is + # Use the year of the last timestep we care about roms_year = time[end_t].year if month == 11: + # For December (note 0-indexed months) best to use the year of the + # first timestep we care about, otherwise we'll have looped to January + # the following year roms_year = time[start_t].year + # Check for leap years + leap_year = False if mod(roms_year, 4) == 0: leap_year = True if mod(roms_year, 100) == 0: @@ -65,10 +100,13 @@ def evap_era_monthly (roms_file, month, colour_bounds=None, save=False, fig_name if leap_year: end_day[1] = 29 - roms_evap = ma.empty(shape(roms_lon)) - roms_evap[:,:] = 0.0 + # Initialise monthly average of ROMS-CICE evaporation + romscice_evap = ma.empty(shape(roms_lon)) + romscice_evap[:,:] = 0.0 num_days = 0 + # Figure out how many days in the 5-day average represented by start_t + # are actually in the month we care about if time[start_t].month-1 == month and time[start_t].day == start_day[month]+2: start_days = 5 elif time[start_t].month-1 == month and time[start_t].day == start_day[month]+1: @@ -83,13 +121,27 @@ def evap_era_monthly (roms_file, month, colour_bounds=None, save=False, fig_name print 'Error: starting index is month ' + str(time[start_t].month) + ', day ' + str(time[start_t].day) return - roms_evap += id.variables['evaporation'][start_t,:-15,:-1]*start_days + # Read ROMS and CICE evaporation, being sure to throw away the ghost cells + # in ROMS, and the northern boundary in both + roms_evap = id_roms.variables['evaporation'][start_t,1:-16,1:-1] + cice_evap = id_cice.variables['evap_ai'][start_t,:-15,:]*cice_conv + # Also read sea ice concentration for scaling + aice = id_cice.variables['aice'][start_t,:-15,:] + # Merge the two evaporation fields based on sea ice concentration + romscice_evap += (aice*cice_evap + (1-aice)*roms_evap)*start_days num_days += start_days + # Between start_t and end_t, we care about all the days for t in range(start_t+1, end_t): - roms_evap += id.variables['evaporation'][t,:-15,:-1]*5 + # As before + roms_evap = id_roms.variables['evaporation'][t,1:-16,1:-1] + cice_evap = id_cice.variables['evap_ai'][t,:-15,:]*cice_conv + aice = id_cice.variables['aice'][t,:-15,:] + romscice_evap += (aice*cice_evap + (1-aice)*roms_evap)*5 num_days += 5 + # Figure out how many days in the 5-day average represented by end_t + # are actually in the month we care about if time[end_t].month-1 == next_month and time[end_t].day == start_day[next_month]+1: end_days = 1 elif time[end_t].month-1 == next_month and time[end_t].day == start_day[next_month]: @@ -104,64 +156,88 @@ def evap_era_monthly (roms_file, month, colour_bounds=None, save=False, fig_name print 'Error: ending index is month ' + str(time[end_t].month) + ', day ' + str(time[end_t].day) return - roms_evap += id.variables['evaporation'][end_t,:-15,:-1]*end_days + # As before + roms_evap = id_roms.variables['evaporation'][end_t,1:-16,1:-1] + cice_evap = id_cice.variables['evap_ai'][end_t,:-15,:]*cice_conv + aice = id_cice.variables['aice'][end_t,:-15,:] + romscice_evap += (aice*cice_evap + (1-aice)*roms_evap)*end_days num_days += end_days + # Make sure we integrated the right number of days if num_days != end_day[month]: print 'Error: found ' + str(num_days) + ' days instead of ' + str(end_day[month]) return - id.close() - roms_evap /= num_days - roms_evap *= 1e6 + id_roms.close() + id_cice.close() + # Convert from integral to average + romscice_evap /= num_days + # Multiply by 1e6 for readability + romscice_evap *= 1e6 + # Read ERA-Interim data for this month and year id = Dataset(era_head + str(roms_year) + era_tail, 'r') + # 1D grid variables era_lon_1d = id.variables['longitude'][:] era_lat_1d = id.variables['latitude'][:] + # Make a 2D version of each era_lon, era_lat = meshgrid(era_lon_1d, era_lat_1d) + # Evaporation data era_evap = id.variables['e'][month,:,:]*era_conv id.close() era_evap *= 1e6 + # 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) era_x = -(era_lat+90)*cos(era_lon*deg2rad+pi/2) era_y = (era_lat+90)*sin(era_lon*deg2rad+pi/2) if colour_bounds is not None: + # User-defined colour bounds lev = linspace(colour_bounds[0], colour_bounds[1], num=50) else: + # Automatic colour bounds min_bound = min(amin(roms_evap), amin(era_evap)) max_bound = max(amax(roms_evap), amax(era_evap)) lev = linspace(min_bound, max_bound, num=50) + # Northern boundary 30S bdry = -30+90 + # Make the figure fig = figure(figsize=(20,9)) + # ERA-Interim fig.add_subplot(1,2,1, aspect='equal') contourf(era_x, era_y, era_evap, lev, extend='both') title('ERA-Interim', fontsize=24) xlim([-bdry, bdry]) ylim([-bdry, bdry]) axis('off') + # MetROMS fig.add_subplot(1,2,2, aspect='equal') - img = contourf(roms_x, roms_y, roms_evap, lev, extend='both') - title('ROMS', fontsize=24) + img = contourf(roms_x, roms_y, romscice_evap, lev, extend='both') + title('MetROMS', fontsize=24) xlim([-bdry, bdry]) ylim([-bdry, bdry]) axis('off') + # Add a colourbar on the bottom cbaxes = fig.add_axes([0.3, 0.04, 0.4, 0.04]) cbar = colorbar(img, orientation='horizontal', cax=cbaxes) + # Main title suptitle(r'Evaporation (10$^{-6}$ kg/m$^2$/s), ' + month_name[month] + ' ' + str(roms_year), fontsize=30) + # Finished if save: fig.savefig(fig_name) else: fig.show() +# Command-line interface if __name__ == "__main__": roms_file = raw_input("Path to ROMS file: ") + cice_file = raw_input("Path to CICE file with the same timesteps as ROMS file: ") month = int(raw_input("Month number (1-12): "))-1 colour_bounds = None get_bounds = raw_input("Set bounds on colour scale (y/n)? ") @@ -176,7 +252,8 @@ def evap_era_monthly (roms_file, month, colour_bounds=None, save=False, fig_name elif action == 'd': save = False fig_name = None - evap_era_monthly(roms_file, month, colour_bounds, save, fig_name) + # Make the plot + evap_era_monthly(roms_file, cice_file, month, colour_bounds, save, fig_name) while True: # Repeat until the user wants to exit @@ -184,7 +261,7 @@ def evap_era_monthly (roms_file, month, colour_bounds=None, save=False, fig_name if repeat == 'y': while True: # Ask for changes to the parameters until the user is done - changes = raw_input("Enter a parameter to change: (1) file path, (2) month number, (3) colour bounds, (4) save/display; or enter to continue: ") + changes = raw_input("Enter a parameter to change: (1) file paths, (2) month number, (3) colour bounds, (4) save/display; or enter to continue: ") if len(changes) == 0: # No more changes to parameters break @@ -192,6 +269,7 @@ def evap_era_monthly (roms_file, month, colour_bounds=None, save=False, fig_name if int(changes) == 1: # New ROMS file roms_file = raw_input("Path to ROMS file: ") + cice_file = raw_input("Path to CICE file with the same timesteps as ROMS file: ") elif int(changes) == 2: # New month month = int(raw_input("Month number (1-12): "))-1 @@ -209,6 +287,6 @@ def evap_era_monthly (roms_file, month, colour_bounds=None, save=False, fig_name # Get a new figure name fig_name = raw_input("File name for figure: ") # Make the plot - evap_era_monthly(roms_file, month, colour_bounds, save, fig_name) + evap_era_monthly(roms_file, cice_file, month, colour_bounds, save, fig_name) else: break diff --git a/file_guide.txt b/file_guide.txt index 68dbf34..b811b23 100644 --- a/file_guide.txt +++ b/file_guide.txt @@ -62,7 +62,7 @@ convert_era.job: A self-submitting batch job which converts 1 year of qsub -v YEAR=1992,COUNT=0 convert_era.job You can also submit multiple years quickly with a bash loop: - for i in `seq 1992 2009`; + for i in `seq 1992 2005`; do qsub -v YEAR=$i,COUNT=0 convert_era.job done @@ -85,7 +85,7 @@ fill_clouds.py: The ERA-Interim dataset for 6-hour total cloud cover has some import *", then "fix_file('filename')" where filename is the path to the ERA-Interim file that has missing timesteps. You can also easily call it in a loop, e.g. - for year in range(1992, 2009+1): + for year in range(1992, 2005+1): fix_file('../data/ERA_Interim/AN_'+str(year)+'_unlim.nc') romscice_atm_smooth.py: At the beginning of spinups, when the ocean is @@ -164,6 +164,18 @@ nbdry_grid_hack.py: Modify the ROMS grid so that the northernmost 15 rows To run: Open python or ipython and type "run nbdry_grid_hack.py". +romscice_tides.py: Create a ROMS tide file containing the first 10 tidal + components interpolated from TPXO 7.2. (This is NOT the same + as Kate's potential tides option; it is for SSH_TIDES rather + than POT_TIDES.) + To run: Make sure the paths near the top of the file, to the + ROMS grid, TPXO 7.2 file, and desired output file, + are correct. Then make sure that you have scipy + version 0.14 or higher (on raijin, this means + switching to python/2.7.6; instructions at the top + of the file). Then open python or ipython and type + "run romscice_tides.py". + add_tide_period.py: Given a ROMS tide file with tidal potential amplitude and phase angle (created using Kate's potential tide scripts), add the tidal period in seconds for each component. @@ -196,8 +208,62 @@ iceberg_melt.py: Read Martin and Adcroft's monthly climatology of freshwater "run iceberg_melt.py". The script will prompt you for the filename of the output ROMS forcing file. - -***POST-PROCESSING DIAGNOSTICS*** +cice_ini.py: Make a CICE initial conditions file (non-standard option + ice_ic='mask') based on NSIDC satellite observations of sea ice + concentration. Wherever concentration exceeds 15%, set it to 100%; + everywhere else, set it to 0%. The CICE code will read this "mask" + and set a constant sea ice thickness of 1 m and snow thickness of + 0.2 m. + To run: Make sure you have an NSIDC file containing the CDR monthly + average of sea ice concentration for the month you want. + Then open python or ipython and type "run cice_ini.py". + The script will prompt you for the path to the NSDIC file, + the ROMS grid file, and the desired output file. + +paul_holland_hack.py: Paul Holland prevents spurious deep convection by applying + an extra 1500 Gt/y of freshwater flux spread evenly south + of 60S. Using the ROMS grid, calculate what the rate of + input should be in kg/m^2/s. Print the value to the + screen. + To run: Open python or ipython and type + "run paul_holland_hack.py". The script will prompt + you for the path to the ROMS grid file. + +romscice_evap.py: Read one year of ERA-Interim evaporation data (12-hourly), + interpolate to the ROMS grid, and add to the existing FC + forcing file created using romscice_atm_subdaily.nc. + To run: Edit user parameters near the top of the script + (mainly just file paths). This script is designed to + be called by a self-submitting batch job, e.g. + era_evap.job. + +era_evap.job: A self-submitting batch job which converts 1 year of ERA-Interim + evaporation 12-hourly data into an existing ROMS-CICE forcing file + (created using romscice_atm_subdaily.py) in chunks of 50 12-hour + timesteps at once (otherwise the python memory overflows because + this is a LOT of data). Depends on romscice_evap.py. + To run: First edit user parameters in romscice_evap.py (see + above). Then edit the PBS job settings at the top of this + file. Then, to get the batch jobs started for a given + year (say 1992, type + qsub -v YEAR=1992,COUNT=0 era_evap.job + You can also submit multiple years quickly with a bash + loop: + for i in `seq 1992 2005`; + do + qsub -v YEAR=$i,COUNT=0 era_evap.job + done + +romscice_gpcp.py: Convert one year of GPCP monthly averages for total + precipitation to a ROMS-CICE input forcing file. + To run: First make sure the paths near the top of the script + (for the ROMS grid, GPCP input files, and desired + output files) are correct. Then open python or ipython + and type "run romscice_gpcp.py". The years 1992-2005 + will be processed in a loop. + + +***POST-PROCESSING DIAGNOSTIC FILES*** make_density_file.py: Given an ocean history or averages file with temperature and salinity data, calculate density fields at each @@ -209,6 +275,21 @@ make_density_file.py: Given an ocean history or averages file with temperature the ocean history or averages file, and the desired path to the new density file. +common_grid.py: Interpolates ROMS output to a regular quarter-degree grid for + easy comparison with FESOM. Writes monthly averages of surface + heat and salt flux, the surface stress vector and its curl, and + the sea ice velocity vector. + To run: First concatenate your ROMS and CICE output into one + long file for each model, containing 5-day averages for + the entire simulation starting from 1 Jan. Then open + python or ipython and type "run common_grid.py". The + script will prompt you for the paths to the ROMS and + CICE output files, and the desired output file on the + common grid. + + +***DIAGNOSTIC FIGURES*** + plot_kinetic_energy.py: Extracts the kinetic energy values written in the ocean.log files (you can supply as many ocean.log files as you like, from e.g. a long run split into several @@ -304,8 +385,16 @@ timeseries_3D.py: Calculates and plots timeseries of total ocean heat content, for the paths to the ROMS grid file, the ocean history or averages file, and the log file. - -***NICE FIGURES*** +aice_animation.py: Create an animation of sea ice concentration for the given + simulation, and save as an mp4 file. + To run: If you are on raijin, first type "module load ffmpeg" + to make sure you will be able to write the mp4 file. + Then edit the variables "directory", "num_ts", and + "start_file" near the top of the file to suit your + simulation. Then open python or ipython and type + "run aice_animation.py". Note that this isn't an + encapsulated function but rather just a script, so + be careful with existing variable names. circumpolar_plot.py: Generates a circumpolar Antarctic plot of the given variable from ROMS. If the variable is depth-dependent, @@ -645,6 +734,21 @@ ssflux_tamura_seasonal.py: Make a figure comparing seasonally-averaged sea ice (and if so, what filename) or display it on the screen. This script is non-repeating. +ssflux_tamura_annual.py: Make a figure comparing annually-averaged sea ice to + ocean salt fluxes, from Tamura's dataset to CICE + output. + To run: First make sure the paths to monthly-averaged, + unit-corrected (to kg/m^2/s) Tamura data files + are correct. Then make sure you have an + annually averaged CICE file. Then open python + or ipython and type + "run ssflux_tamura_annual.py". The script will + prompt you for the path to this CICE file, the + year it corresponds to, and whether you want + to save the figure (and if so, what filename) + or display it on the screen. This script + is non-repeating. + sose_seasonal.py: Make a 4x2 plot showing lat vs. depth slices of seasonally averaged temperature (top) and salinity (bottom) at the given longitude, for the 2005-2010 SOSE climatology. @@ -702,38 +806,83 @@ seaice_budget_thermo.py: Create four plots showing timeseries of the so, what filenames) or display them on the screen. - -***ANIMATIONS*** - -aice_animation.py: Create an animation of sea ice concentration for the given - simulation, and save as an mp4 file. - To run: If you are on raijin, first type "module load ffmpeg" - to make sure you will be able to write the mp4 file. - Then edit the variables "directory", "num_ts", and - "start_file" near the top of the file to suit your - simulation. Then open python or ipython and type - "run aice_animation.py". Note that this isn't an - encapsulated function but rather just a script, so - be careful with existing variable names. +compare_nic_seasonal.py: Make a 4x2 plot comparing seasonal averages of the + given variable from my implementation of CICE in + MetROMS (bottom) and Nic Hannah's implementation + coupled to MOM (top). + To run: Open python or ipython and type + "run compare_nic_seasonal.py". The script will + prompt you for the path to the CICE-MetROMS + output file, the variable name to plot, + optional colour bounds, and whether you want + to save the figure (and if so, what filename) + or display it on the screen. The script will + repeat as many times as you like. + +era_gpcp.py: Plot precipitation (1992-2005 average) in the ERA-Interim versus + GPCP datasets over the Southern Ocean. + To run: First make sure the ERA-Interim and GPCP paths near the top + of the file are correct. Then open python or ipython and + type "run era_gpcp.py". The figure will be saved under + the filename "era_gpcp.png". + +evap_era_monthly.py: Plot monthly-averaged evaporation as calculated by MetROMS + (ROMS and CICE calculate separately, so merge them based + on sea ice concentration) and compare it to ERA-Interim. + To run: First make sure the ERA-Interim paths near the top + of the file are correct. Then open python or + ipython and type "run evap_era_monthly.py". The + script will prompt you for paths to the ROMS and + CICE output files (5-day averages for the same + timesteps; if multiple instances of the given month + are present the script will select the last one), + the month to plot, optional colour bounds, and + whether you want to save the figure (and if so, + what filename) or display it on the screen. The + script will repeat as many times as you like. + +holland_fig1.py: Recreate Figure 1 of Holland et al 2014 + (doi:10.1175/JCLI-D-13-00301.1) using ROMS output: annually + averaged barotropic streamfunction, JJA mean mixed-layer depth + (just plot the KPP boundary layer depth for now), annually + averaged bottom water temperature and salinity. + To run: Open python or ipython and type "run holland_fig1.py". + The script will prompt you for the path to the ROMS + grid file, and the ROMS output file containing exactly + 1 year of 5-day averages including a continuous + June-August period. The figure will be displayed on + the screen. + +temp_salt_slice.py: Create a 2x1 plot showing zonal slices (depth vs latitude) + of temperature and salinity interpolated to the given + longitude, at the given timestep. + To run: Open python or ipython and type + "run temp_salt_slice.py". The script will prompt you + for the path to the ROMS averages file, the timestep + to plot, the longitude to interpolate to, the + deepest depth to plot, and whether you want to save + the figure (and if so, what filename) or display it + on the screen. This script repeats as many times + as you want. ***FIGURES FOR ADVECTION PAPER*** -adv_amery_tsplots_indiv.py: For each advection experiment, plot zonal slices of - temperature and salinity through 71E (Amery Ice - Shelf) at the end of the simulation. This is really - just a special case of zonal_plot.py. - To run: Make sure the paths to each simulation are - correct. Then open python or ipython and - type "run adv_amery_tsplots_indiv.py". - adv_amery_tsplots.py: Create a 2x2 plot showing zonal slices of temperature and salinity through 71E (Amery Ice Shelf) at the end of the U3_LIM and C4_LD simulations. To run: Make sure the paths to the simulation directories and 31 December ROMS files are correct. Then open python or ipython and type - "run adv_amery_tsplots.py". + "run adv_amery_tsplots.py". The figure will be + saved under the filename adv_amery_tsplots.png. + +adv_ross_tsplots.py: Same as adv_amery_tsplots.py, but for the Ross Sea (180E). + To run: Make sure the paths to the simulation directories + and 31 December ROMS files are correct. Then open + python or ipython and type + "run adv_ross_tsplots.py". The figure will be + saved under the filename adv_ross_tsplots.png. adv_freezingpt_slice.py: Plot the difference from the freezing temperature at a specific timestep through the Weddell Sea in the @@ -743,7 +892,9 @@ adv_freezingpt_slice.py: Plot the difference from the freezing temperature at really just a special case of freezingpt_slice.py. To run: Make sure the file_path is correct. Then open python or ipython and type - "run adv_freezingpt_slice.py". + "run adv_freezingpt_slice.py". The figure will + be saved under the filename + adv_freezingpt_slice.png. adv_timeseries_volume.py: Plot timeseries of total sea ice volume for all the advection experiments. Before running this script, @@ -751,39 +902,109 @@ adv_timeseries_volume.py: Plot timeseries of total sea ice volume for all the To run: Make sure the paths to directories and the sea ice logfiles are correct. Then open python or ipython and type - "run adv_timeseries_volume.py". + "run adv_timeseries_volume.py". The figure + will be saved under the filename + adv_timeseries_volume.png. adv_frazil.py: Create a 3x2 plot of annually averaged frazil formation for each advection experiment: the absolute values for U3_LIM, and the anomalies from U3_LIM for the other 5 experiments. To run: Make sure the paths to the simulation directories and annually averaged CICE files are correct. Then open - python or ipython and type "run adv_frazil.py". + python or ipython and type "run adv_frazil.py". The + figure will be saved under the filename adv_frazil.png. adv_thickness.py: Like adv_frazil.py, but for sea ice effective thickness on 23 August (the sea ice area max). To run: Make sure the paths to the simulation directories and 23 August CICE files are correct. Then open python or - ipython and type "run adv_frazil.py". + ipython and type "run adv_thickness.py". The figure + will be saved under the filename adv_thickness.png. adv_sss.py: Like adv_thickness.py, but for sea surface salinity. To run: Make sure the paths to the simulation directories and 23 August CICE files are correct. Then open python or ipython - and type "run adv_sss.py". + and type "run adv_sss.py". The figure will be saved under + the filename adv_sss.png. adv_mld.py: Like adv_thickness.py, but for mixed layer depth as calculated by the KPP parameterisation. To run: Make sure the paths to the simulation directories and 23 August ROMS files are correct. Then open python or ipython - and type "run adv_mld.py". + and type "run adv_mld.py". The figure will be saved under + the filename adv_mld.png. -adv_polynyas.py: Create a 1x2 plot of sea ice concentration along the coast of - Wilkes Land on 23 August (the sea ice area max), showing the +adv_polynyas.py: Create a 1x2 plot of sea ice concentration near the Amery Ice + Shelf on 23 August (the sea ice area max), showing the difference in polynya size between the U3_LIM and C4_LD simulations. To run: Make sure the paths to the simulation directories and 23 August CICE files are correct. Then open python or - ipython and type "run adv_polynyas.py". + ipython and type "run adv_polynyas.py". The figure + will be saved under the filename adv_polynyas.png. + +adv_u3_sss_anom.py: Plot the sea surface salinity anomaly between the U3 and + U3_LIM simulations on 23 August (sea ice area maximum). + To run: Make sure the paths to the simulation directories + and 23 August CICE files are correct. Then open + python or ipython and type "run adv_u3_sss_anom.py". + The figure will be saved under the filename + sss_u3_anom.png. + +timeseries_sss.py: Calculates and plots timeseries of area-averaged sea surface + salinity and surface salt flux during a ROMS simulation. + Also writes the timeseries to a log file so it doesn't have + to be recomputed if the run is extended; at the beginning of + this script, previous values will be read from the same log + file if it exists. + To run: Open python or ipython and type + "run timeseries_sss.py". The script will prompt you + for the path to the ROMS averages file and the log + file. + +timeseries_i2osalt.py: Calculates and plots timeseries of area-averaged sea ice + to ocean salt flux during a ROMS simulation. Also writes + the timeseries to a log file so it doesn't have to be + recomputed if the run is extended; at the beginning of + this script, previous values will be read from the same + log file if it exists. + To run: Open python or ipython and type + "run timeseries_i2osalt.py". The script will + prompt you for the path to the CICE history file + and the log file. + + +***FIGURES FOR GARY'S LIMITERS PAPER*** + +gadv_freezingpt_slice.py: As adv_freezingpt_slice.py, but for the U3 experiment + through a different slice (through 8E) and timestep. + To run: Make sure the path to the ocean history file + is correct. Then open python or ipython and + type "run gadv_freezingpt_slice.py". The + figure will be saved with the filename + "kn_fig1.png". + +gadv_frazil_thickness.py: Plot anomalies between the U3 and U3_LIM simulations + in 2 circumpolar plots. On the left, annually + averaged frazil formation (as in adv_frazil.py). On + the right, effective sea ice thickness on 23 August + (sea ice area max). + To run: Make sure the paths to the simulation + directories near the top of the file are + correct. Then open python or ipython and type + "run gadv_frazil_thickness.py". The figure + will be saved with the filename "kn_fig2.png". + +gadv_frazil_bathy.py: Plot a section of the coastline from 25-45E. On the left, + show anomalies in annually averaged frazil formation + between the U3 and U3_LIM simulations. On the right, show + bathymetry. This way the relationship between spurious + supercooling and steep bathymetry is obvious. + To run: Make sure the paths to the simulation directories + and ROMS grid file near the top of the file are + correct. Then open python or ipython and type + "run gadv_frazil_bathy.py". The figure will be + saved with the filename "kn_fig3.png". ***UTILITY FUNCTIONS*** @@ -824,6 +1045,49 @@ unesco.py: Calculates the UNESCO seawater equation of state (1980): given To run: The function unesco is designed to be called by another script. See make_density_file.py for an example. +interp_era2roms.py: Given an array on the ERA-Interim grid, interpolate any + missing values, and then interpolate to the ROMS grid. + To run: The function interp_era2roms is designed to be + called by another script. See + romscice_atm_subdaily.py for an example. + +interp_lon_roms.py: Linearly interpolate ROMS data, z, and latitude to the + specified longitude. + To run: The function interp_lon_roms is designed to be + called by another script. See zonal_plot.py for + an example. + +interp_lon_sose.py: Linearly interpolate SOSE data to the specified longitude. + To run: The function interp_lon_sose is designed to be + called by another script. See sose_roms_seasonal.py + for an example. + +monthly_avg_cice.py: Average the given variable in the given CICE file over the + given month. + To run: The function monthly_avg_cice is designed to be + called by another script. See nsidc_aice_monthly.py + for an example. + +monthly_avg_roms.py: Average the given variable in the given ROMS file over the + given month. + To run: The function monthly_avg_roms is designed to be + called by another script. See common_grid.py + for an example. + +seasonal_avg_cice.py: Calculate seasonal averages (DJF, MAM, JJA, SON) of the + given variable in the given CICE file. + To run: The function seasonal_avg_cice is designed to be + called by another script. See aice_hi_seasonal.py + for an example. + +seasonal_avg_roms.py: Calculate seasonal averages (DJF, MAM, JJA, SON) of the + given variable in the given ROMS file. + To run: The function seasonal_avg_roms is designed to be + called by another script. See + temp_salt_seasonal.py for an example. + + + diff --git a/gadv_frazil_bathy.py b/gadv_frazil_bathy.py index 6b24c65..1dd0002 100644 --- a/gadv_frazil_bathy.py +++ b/gadv_frazil_bathy.py @@ -2,49 +2,67 @@ from numpy import * from matplotlib.pyplot import * +# Plot a section of the coastline from 25-45E. On the left, show anomalies +# in annually averaged frazil formation between the U3 and U3_LIM simulations. +# On the right, show bathymetry. This way the relationship between spurious +# supercooling and steep bathymetry is obvious. def gadv_frazil_bathy (): + # Paths to simulation directories path_up3l = '/short/m68/kaa561/advection_bitreproducible/u3_lim/' path_up3 = '/short/m68/kaa561/advection_bitreproducible/u3/' + # Path to ROMS grid (containing bathymetry data) roms_grid = '/short/m68/kaa561/metroms_iceshelf/apps/common/grid/circ30S_quarterdegree.nc' + # Filename for annually averaged CICE data file_tail = 'cice/rundir/history/iceh_avg.nc' + # Bounds on lat and lon lon_min = 25 lon_max = 45 lat_min = -71 lat_max = -64 + # Bounds for colour scales max_frazil = 0.4 tick_frazil = 0.2 max_bathy = 5 tick_bathy = 1 + # Read U3_LIM grid and frazil data id = Dataset(path_up3l + file_tail, 'r') lon_cice = id.variables['TLON'][50:250,0:200] lat_cice = id.variables['TLAT'][50:250,0:200] frazil0 = id.variables['frazil'][0,50:250,0:200] id.close() + # Read U3 frazil data id = Dataset(path_up3 + file_tail, 'r') frazil1 = id.variables['frazil'][0,50:250,0:200] id.close() + # Calculate anomaly frazil_anom = frazil1 - frazil0 + # Read bathymetry id = Dataset(roms_grid, 'r') lon_roms = id.variables['lon_rho'][51:251,1:201] lat_roms = id.variables['lat_rho'][51:251,1:201] bathy = id.variables['h'][51:251,:201]*1e-3 + # Mask out land and ice shelf points mask = id.variables['mask_rho'][51:251,:201]-id.variables['mask_zice'][51:251,:201] id.close() bathy = ma.masked_where(mask==0, bathy) + # Make the plot fig = figure(figsize=(18,8)) + # Frazil anomalies ax = fig.add_subplot(1, 2, 1) img1 = pcolor(lon_cice, lat_cice, frazil_anom, vmin=-max_frazil, vmax=max_frazil, cmap='RdBu_r') title('a) Frazil ice formation (cm/day), annual average\nAnomalies: UP3 minus UP3L', fontsize=20) xlabel('Longitude', fontsize=18) xlim([lon_min, lon_max]) ylim([lat_min, lat_max]) + # Colour bar cbaxes1 = fig.add_axes([0.05, 0.25, 0.015, 0.5]) cbar1 = colorbar(img1, ticks=arange(-max_frazil, max_frazil+tick_frazil, tick_frazil), cax=cbaxes1, extend='both') cbar1.ax.tick_params(labelsize=16) + # Make sure lon and lat ticks are the way we want them lon_ticks = arange(lon_min, lon_max+5, 5) ax.set_xticks(lon_ticks) lon_labels = [] @@ -55,8 +73,10 @@ def gadv_frazil_bathy (): ax.set_yticks(lat_ticks) lat_labels = [] for val in lat_ticks: + # No latitude labels on this plot, they bump into the colourbar lat_labels.append('') ax.set_yticklabels(lat_labels, fontsize=16) + # Bathymetry ax = fig.add_subplot(1, 2, 2) img2 = pcolor(lon_roms, lat_roms, bathy, vmin=0, vmax=max_bathy, cmap='jet') title('b) Bathymetry (km)', fontsize=20) @@ -84,6 +104,7 @@ def gadv_frazil_bathy (): fig.savefig('kn_fig3.png') +# Command-line interface if __name__ == "__main__": gadv_frazil_bathy() diff --git a/gadv_frazil_thickness.py b/gadv_frazil_thickness.py index 32f7e95..b3d3f78 100644 --- a/gadv_frazil_thickness.py +++ b/gadv_frazil_thickness.py @@ -2,13 +2,20 @@ from numpy import * from matplotlib.pyplot import * +# Plot anomalies between the U3 and U3_LIM simulations in 2 circumpolar plots. +# On the left, annually averaged frazil formation. On the right, effective sea +# ice thickness on 23 August (sea ice area max). def gadv_frazil_thickness (): + # Paths to simulation directories path_up3l = '/short/m68/kaa561/advection_bitreproducible/u3_lim/' path_up3 = '/short/m68/kaa561/advection_bitreproducible/u3/' + # Annually averaged CICE file file_tail_avg = 'cice/rundir/history/iceh_avg.nc' + # Daily averaged CICE file for 23 August file_tail_23aug = 'cice/rundir/history/iceh.1992-08-23.nc' + # Bounds for colour scales max_frazil = 0.5 tick_frazil = 0.25 max_thickness = 1 @@ -24,11 +31,16 @@ def gadv_frazil_thickness (): # Boundary of regular grid to embed circle in circle_bdry = -70+90 + # Labels for longitude around the plot + # Lat and lon locations of labels lon_ticks = array([-120, -60, 60, 120, 180]) lat_ticks = array([-58, -56, -54, -55.5, -56]) + # Label text lon_labels = [r'120$^{\circ}$W', r'60$^{\circ}$W', r'60$^{\circ}$E', r'120$^{\circ}$E', r'180$^{\circ}$'] + # Rotation of text lon_rot = [-60, 60, -60, 60, 0] + # Read U3_LIM frazil data and grid id = Dataset(path_up3l + file_tail_avg, 'r') frazil_tmp = id.variables['frazil'][0,:275,:] lon_tmp = id.variables['TLON'][:275,:] @@ -49,21 +61,23 @@ def gadv_frazil_thickness (): frazil0[:,:-1] = frazil_tmp frazil0[:,-1] = frazil_tmp[:,0] + # Read U3 frazil data id = Dataset(path_up3 + file_tail_avg, 'r') frazil_tmp = id.variables['frazil'][0,:275,:] id.close() frazil1 = ma.empty([size(frazil_tmp,0), size(frazil_tmp,1)+1]) frazil1[:,:-1] = frazil_tmp frazil1[:,-1] = frazil_tmp[:,0] + # Calculate anomaly frazil_anom = frazil1 - frazil0 + # Repeat for thickness data id = Dataset(path_up3l + file_tail_23aug, 'r') thickness_tmp = id.variables['aice'][0,:275,:]*id.variables['hi'][0,:275,:] id.close() thickness0 = ma.empty([size(thickness_tmp,0), size(thickness_tmp,1)+1]) thickness0[:,:-1] = thickness_tmp thickness0[:,-1] = thickness_tmp[:,0] - id = Dataset(path_up3 + file_tail_23aug, 'r') thickness_tmp = id.variables['aice'][0,:275,:]*id.variables['hi'][0,:275,:] id.close() @@ -90,7 +104,9 @@ def gadv_frazil_thickness (): 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) + # Make the plot fig = figure(figsize=(18,8)) + # Frazil ax = fig.add_subplot(1, 2, 1, aspect='equal') # First shade land contourf(x, y, land, 1, colors=(('0.6', '0.6', '0.6'))) @@ -102,10 +118,11 @@ def gadv_frazil_thickness (): text(x_ticks[i], y_ticks[i], lon_labels[i], ha='center', rotation=lon_rot[i], fontsize=15) axis('off') title('a) Frazil ice formation (cm/day), annual average', fontsize=20) + # Colourbar cbaxes1 = fig.add_axes([0.07, 0.25, 0.015, 0.5]) cbar1 = colorbar(img1, ticks=arange(-max_frazil, max_frazil+tick_frazil, tick_frazil), cax=cbaxes1, extend='both') cbar1.ax.tick_params(labelsize=16) - + # Thickness ax = fig.add_subplot(1, 2, 2, aspect='equal') contourf(x, y, land, 1, colors=(('0.6', '0.6', '0.6'))) contourf(x_reg, y_reg, land_circle, 1, colors=(('0.6', '0.6', '0.6'))) @@ -118,6 +135,7 @@ def gadv_frazil_thickness (): cbar2 = colorbar(img2, ticks=arange(-max_thickness, max_thickness+tick_thickness, tick_thickness), cax=cbaxes2, extend='both') cbar2.ax.tick_params(labelsize=16) + # Main title suptitle('Anomalies: UP3 minus UP3L', fontsize=26) subplots_adjust(wspace=0.05) @@ -125,6 +143,7 @@ def gadv_frazil_thickness (): fig.savefig('kn_fig2.png') +# Command-line interface if __name__ == "__main__": gadv_frazil_thickness() diff --git a/gadv_freezingpt_slice.py b/gadv_freezingpt_slice.py index b64b07c..f8747b7 100644 --- a/gadv_freezingpt_slice.py +++ b/gadv_freezingpt_slice.py @@ -3,6 +3,9 @@ from matplotlib.pyplot import * from calc_z import * +# Plot the difference from the freezing temperature at a specific timestep +# through an i-slice near 8E in the U3 experiment. There is no time-averaging, +# spatial averaging, or interpolation. This shows off the spurious supercooling. def gadv_freezingpt_slice (): # Path to ocean history file @@ -62,6 +65,7 @@ def gadv_freezingpt_slice (): xlim([lat_min, lat_max]) ylim([depth_min, 0]) + # Make sure lat and lon ticks are the way we want them lat_ticks = arange(lat_min+0.5, lat_max+0.5, 1) xticks(lat_ticks) lat_labels = [] diff --git a/holland_fig1.py b/holland_fig1.py index 46c2d63..bd2cf8e 100644 --- a/holland_fig1.py +++ b/holland_fig1.py @@ -4,6 +4,14 @@ from cartesian_grid_2d import * from rotate_vector_roms import * +# Recreate Figure 1 of Holland et al 2014 (doi:10.1175/JCLI-D-13-00301.1) using +# ROMS output: annually averaged barotropic streamfunction, JJA mean mixed-layer +# depth (just plot the KPP boundary layer depth for now), annually averaged +# bottom water temperature and salinity. +# Input: +# grid_path = path to ROMS grid file +# file_path = path to ROMS output file containing 1 year of 5-day averages, +# including a continuous June-August period def holland_fig1 (grid_path, file_path): deg2rad = pi/180.0 @@ -23,19 +31,29 @@ def holland_fig1 (grid_path, file_path): fig = figure(figsize=(16,12)) # Barotropic streamfunction + # First read bartropic velocity vector id = Dataset(file_path, 'r') ubar_xy = mean(id.variables['ubar'][:,:-15,:], axis=0) vbar_xy = mean(id.variables['vbar'][:,:-15,:], axis=0) id.close() + # Rotate to lon-lat space ubar, vbar = rotate_vector_roms(ubar_xy, vbar_xy, angle) + # Throw away the overlapping periodic boundary ubar = ubar[:,1:] + # Mask ice shelves ubar = ma.masked_where(zice!=0, ubar) + # Water column thickness wct = h+zice + # Horizontal differentials dx, dy = cartesian_grid_2d(lon, lat) + # Indefinite integral from south to north of u*dz*dy, convert to Sv baro_strf = cumsum(ubar*wct*dy, axis=0)*1e-6 + # Colour levels lev1 = arange(-50, 150+10, 10) + # Plot ax1 = fig.add_subplot(2, 2, 1, aspect='equal') img = contourf(x, y, baro_strf, lev1, extend='both') + # Contour 0 Sv in black contour(x, y, baro_strf, levels=[0], colors=('black')) title('Barotropic streamfunction (Sv)', fontsize=24) xlim([-35, 39]) @@ -46,17 +64,19 @@ def holland_fig1 (grid_path, file_path): cbar1.ax.tick_params(labelsize=16) # JJA mixed layer depth - start_month = 6 - end_month = 8 - start_day = 1 - next_startday = 1 - end_day = 31 - prev_endday = 31 - ndays_season = 92 + start_month = 6 # Start in June + end_month = 8 # End in August + start_day = 1 # First day in June + next_startday = 1 # First day in September + end_day = 31 # Last day in August + prev_endday = 31 # Last day in May + ndays_season = 92 # Number of days in June+July+August id = Dataset(file_path, 'r') + # Read time axis and get dates time_id = id.variables['ocean_time'] time = num2date(time_id[:], units=time_id.units, calendar=time_id.calendar.lower()) - end_t = -1 + # Find the last timestep we care about + end_t = -1 # Missing value flag for t in range(size(time)-1, -1, -1): if time[t].month == end_month and time[t].day in range(end_day-2, end_day+1): end_t = t @@ -64,10 +84,12 @@ def holland_fig1 (grid_path, file_path): if time[t].month == end_month+1 and time[t].day in range(next_startday, next_startday+2): end_t = t break + # Make sure we actually found it if end_t == -1: print 'Error: ' + file_path + ' does not contain a complete JJA' return - start_t = -1 + # Find the first timestep we care about + start_t = -1 # Missing value flag for t in range(end_t, -1, -1): if time[t].month == start_month-1 and time[t].day in range(prev_endday-1, prev_endday+1): start_t = t @@ -75,12 +97,15 @@ def holland_fig1 (grid_path, file_path): if time[t].month == start_month and time[t].day in range(start_day, start_day+3): start_t = t break + # Make sure we found it if start_t == -1: print 'Error: ' + file_path + ' does not contain a complete JJA' return + # Initialise time-averaged KPP boundary layer depth hsbl = ma.empty(shape(lon)) - hsbl[:,:] = 0.0 + hsbl[:,:] = 0.0 ndays = 0 + # Figure out how many of the 5 days represented in start_t we care about if time[start_t].month == start_month and time[start_t].day == start_day+2: start_days = 5 elif time[start_t].month == start_month and time[start_t].day == start_day+1: @@ -94,11 +119,14 @@ def holland_fig1 (grid_path, file_path): else: print 'Error: starting index is month ' + str(time[start_t].month) + ', day ' + str(time[start_t].day) return + # Integrate Hsbl weighted by start_days hsbl += id.variables['Hsbl'][start_t,:-15,1:]*start_days ndays += start_days + # Between start_t and end_t, we care about all the days for t in range(start_t+1, end_t): hsbl += id.variables['Hsbl'][t,:-15,1:]*5 ndays += 5 + # Figure out how many of the 5 days represented in end_t we care about if time[end_t].month == end_month+1 and time[end_t].day == next_startday+1: end_days = 1 elif time[end_t].month == end_month+1 and time[end_t].day == next_startday: @@ -112,17 +140,23 @@ def holland_fig1 (grid_path, file_path): else: print 'Error: ending index is month ' + str(time[end_t].month) + ', day ' + str(time[end_t].day) return + # Integrate weighted by end_days hsbl += id.variables['Hsbl'][end_t,:-15,1:]*end_days ndays += end_days if ndays != ndays_season: print 'Error: found ' + str(ndays) + ' days instead of ' + str(ndays_season) return id.close() + # Convert from integral to average hsbl[:,:] = hsbl[:,:]/ndays + # Mask out ice shelves, change sign, and call it mixed layer depth mld = ma.masked_where(zice!=0, -hsbl) + # Colour levels lev2 = arange(0, 300+25, 25) + # Plot ax2 = fig.add_subplot(2, 2, 2, aspect='equal') img = contourf(x, y, mld, lev2, extend='both') + # Contour 100 m in black contour(x, y, mld, levels=[100], colors=('black')) title('Winter mixed layer depth (m)', fontsize=24) xlim([-35, 39]) @@ -136,10 +170,14 @@ def holland_fig1 (grid_path, file_path): id = Dataset(file_path, 'r') bwtemp = mean(id.variables['temp'][:,0,:-15,1:], axis=0) id.close() + # Mask ice shelves bwtemp = ma.masked_where(zice!=0, bwtemp) + # Colour levels lev3 = arange(-2, 2+0.2, 0.2) + # Plot ax3 = fig.add_subplot(2, 2, 3, aspect='equal') img = contourf(x, y, bwtemp, lev3, extend='both') + # Contour 0C in black contour(x, y, bwtemp, levels=[0], colors=('black')) title(r'Bottom temperature ($^{\circ}$C', fontsize=24) xlim([-35, 39]) @@ -157,6 +195,7 @@ def holland_fig1 (grid_path, file_path): lev4 = arange(34.5, 34.8+0.025, 0.025) ax4 = fig.add_subplot(2, 2, 4, aspect='equal') img = contourf(x, y, bwsalt, lev4, extend='both') + # Contour 34.65 psu in black contour(x, y, bwsalt, levels=[34.65], colors=('black')) title('Bottom salinity (psu)', fontsize=24) xlim([-35, 39]) @@ -170,6 +209,7 @@ def holland_fig1 (grid_path, file_path): #fig.savefig('fig1_year3.png') +# Command-line interface if __name__ == "__main__": grid_path = raw_input("Path to ROMS grid file: ") diff --git a/monthly_avg_cice.py b/monthly_avg_cice.py index 014d0f7..4d218a2 100644 --- a/monthly_avg_cice.py +++ b/monthly_avg_cice.py @@ -1,6 +1,18 @@ from netCDF4 import Dataset, num2date from numpy import * +# Average the given variable in the given CICE file over the given month. +# Input: +# file_path = path to CICE output file containing 5-day averages, including at +# least one complete instance of the given month. +# var = variable name +# shape = vector containing the dimensions (excluding time) of the variable +# month = month to average over (0-11) +# instance = optional integer indicating which instance of the given month in +# this file we should use. For instance=1 use the first instance, +# etc. If instance=-1 (the default) the last instance is used. +# Output: +# monthly_data = array of data averaged over the given month. def monthly_avg_cice (file_path, var, shape, month, instance=-1): # Starting and ending days in each month @@ -62,10 +74,15 @@ def monthly_avg_cice (file_path, var, shape, month, instance=-1): print 'Error: ' + file_path + ' does not contain ' + str(instance) + ' ' + month_name[month] + 's' return - leap_year = False + # Figure out what year it is, based on the year of the last timestep we care + # about cice_year = cice_time[end_t].year if month == 11: + # December: the last timestep might be January next year + # Use the year of the first timestep we care about cice_year = cice_time[start_t].year + # Check for leap years + leap_year = False if mod(cice_year, 4) == 0: leap_year = True if mod(cice_year, 100) == 0: @@ -98,9 +115,11 @@ def monthly_avg_cice (file_path, var, shape, month, instance=-1): print 'Error: starting index is month ' + str(cice_time[start_t].month) + ', day ' + str(cice_time[start_t].day) return + # Integrate data weighted by start_days monthly_data += id.variables[var][start_t,:]*start_days num_days += start_days + # Between start_t and end_t, we care about all the days for t in range(start_t+1, end_t): monthly_data += id.variables[var][t,:]*5 num_days +=5 @@ -126,14 +145,17 @@ def monthly_avg_cice (file_path, var, shape, month, instance=-1): print 'Error: ending index is month ' + str(cice_time[end_t].month) + ', day ' + str(cice_time[end_t].day) return + # Integrate data weighted by end_days monthly_data += id.variables[var][end_t,:]*end_days num_days += end_days + # Make sure we got the right number of days if num_days != end_day[month]: print 'Error: found ' + str(num_days) + ' days instead of ' + str(end_day[month]) return id.close() + # Convert from integral to average monthly_data /= num_days return monthly_data diff --git a/monthly_avg_roms.py b/monthly_avg_roms.py index 3944770..3ed5f49 100644 --- a/monthly_avg_roms.py +++ b/monthly_avg_roms.py @@ -1,6 +1,18 @@ from netCDF4 import Dataset, num2date from numpy import * +# Average the given variable in the given ROMS file over the given month. +# Input: +# file_path = path to ROMS output file containing 5-day averages, including at +# least one complete instance of the given month. +# var = variable name +# shape = vector containing the dimensions (excluding time) of the variable +# month = month to average over (0-11) +# instance = optional integer indicating which instance of the given month in +# this file we should use. For instance=1 use the first instance, +# etc. If instance=-1 (the default) the last instance is used. +# Output: +# monthly_data = array of data averaged over the given month. def monthly_avg_roms (file_path, var, shape, month, instance=-1): # Starting and ending days in each month @@ -77,10 +89,15 @@ def monthly_avg_roms (file_path, var, shape, month, instance=-1): print 'Error: ' + file_path + ' does not contain ' + str(instance) + ' ' + month_name[month] + 's' return - leap_year = False + # Figure out what year it is, based on the year of the last timestep we care + # about roms_year = time[end_t].year if month == 11: + # December: the last timestep might be January next year + # Use the year of the first timestep we care about roms_year = time[start_t].year + # Check for leap years + leap_year = False if mod(roms_year, 4) == 0: leap_year = True if mod(roms_year, 100) == 0: @@ -92,6 +109,8 @@ def monthly_avg_roms (file_path, var, shape, month, instance=-1): num_days = 0 + # Figure out how many of the 5 days averaged in start_t are actually within + # this month if time[start_t].month-1 == month and time[start_t].day == start_day[month]+2: start_days = 5 elif time[start_t].month-1 == month and time[start_t].day == start_day[month]+1: @@ -106,13 +125,17 @@ def monthly_avg_roms (file_path, var, shape, month, instance=-1): print 'Error: starting index is month ' + str(time[start_t].month) + ', day ' + str(time[start_t].day) return + # Integrate data weighted by start_days monthly_data += id.variables[var][start_t,:]*start_days num_days += start_days + # Between start_t and end_t, we care about all the days for t in range(start_t+1, end_t): monthly_data += id.variables[var][t,:]*5 num_days += 5 + # Figure out how many of the 5 days averaged in end_t are actually within + # this month if time[end_t].month-1 == next_month and time[end_t].day == start_day[next_month]+1: end_days = 1 elif time[end_t].month-1 == next_month and time[end_t].day == start_day[next_month]: @@ -127,14 +150,17 @@ def monthly_avg_roms (file_path, var, shape, month, instance=-1): print 'Error: ending index is month ' + str(time[end_t].month) + ', day ' + str(time[end_t].day) return + # Integrate data weighted by end_days monthly_data += id.variables[var][end_t,:]*end_days num_days += end_days + # Make sure we got the right number of days if num_days != end_day[month]: print 'Error: found ' + str(num_days) + ' days instead of ' + str(end_day[month]) return id.close() + # Convert from integral to average monthly_data /= num_days return monthly_data diff --git a/nbdry_grid_hack.py b/nbdry_grid_hack.py index e327df7..c8b569e 100644 --- a/nbdry_grid_hack.py +++ b/nbdry_grid_hack.py @@ -42,6 +42,6 @@ def nbdry_grid_hack (grid_file, num_pts): # Command-line interface if __name__ == "__main__": - grid_file='../ROMS-CICE-MCT/apps/common/grid/circ30S_quarterdegree_good.nc' + grid_file='../metroms_iceshelf/apps/common/grid/circ30S_quarterdegree.nc' num_pts = 15 nbdry_grid_hack(grid_file, num_pts) diff --git a/nic_budget_ross_weddell.py b/nic_budget_ross_weddell.py deleted file mode 100644 index 1b67e0d..0000000 --- a/nic_budget_ross_weddell.py +++ /dev/null @@ -1,131 +0,0 @@ -from netCDF4 import Dataset -from numpy import * -from matplotlib.pyplot import * -from matplotlib.font_manager import FontProperties -from cartesian_grid_2d import * - -def nic_budget_ross_weddell (loc_flag, save=False, fig_name=None): - - nic_dir_head = '/g/data/gh5/access_om_025-CORE_NYF/output' - output_number = arange(135, 137+1) - nic_dir_tail = '/ice/HISTORY/' - nic_file_head = 'iceh.0' - nic_year_number = output_number-4 - days_per_month = [31,28,31,30,31,30,31,31,30,31,30,31] - num_time = 12*size(output_number) - - if loc_flag == 'r': - min_i = 300 - max_i = 600 - min_j = 30 - max_j = 120 - elif loc_flag == 'w': - min_i = 860 - max_i = 1100 - min_j = 30 - max_j = 140 - - # Read CICE grid - id = Dataset(nic_dir_head + str(output_number[0]) + nic_dir_tail + nic_file_head + str(nic_year_number[0]) + '-01.nc', 'r') - lon = id.variables['TLON'][:,:] - lat = id.variables['TLAT'][:,:] - id.close() - - # Calculate elements of area - dx, dy = cartesian_grid_2d(lon, lat) - dA = dx*dy - dA = dA[min_j:max_j, min_i:max_i] - - time = arange(num_time)/12.0 - - fontP = FontProperties() - fontP.set_size('small') - - aice = ma.empty([size(time), size(dA,0), size(dA,1)]) - congel = ma.empty([size(time), size(dA,0), size(dA,1)]) - frazil = ma.empty([size(time), size(dA,0), size(dA,1)]) - snoice = ma.empty([size(time), size(dA,0), size(dA,1)]) - meltt = ma.empty([size(time), size(dA,0), size(dA,1)]) - meltb = ma.empty([size(time), size(dA,0), size(dA,1)]) - meltl = ma.empty([size(time), size(dA,0), size(dA,1)]) - - t = 0 - for year in range(size(output_number)): - for month in range(12): - if month+1 < 10: - filename = nic_dir_head + str(output_number[year]) + nic_dir_tail + nic_file_head + str(nic_year_number[year]) + '-0' + str(month+1) + '.nc' - else: - filename = nic_dir_head + str(output_number[year]) + nic_dir_tail + nic_file_head + str(nic_year_number[year]) + '-' + str(month+1) + '.nc' - id = Dataset(filename, 'r') - aice[t,:,:] = id.variables['aice'][0,min_j:max_j,min_i:max_i] - congel[t,:,:] = id.variables['congel'][0,min_j:max_j,min_i:max_i] - frazil[t,:,:] = id.variables['frazil'][0,min_j:max_j,min_i:max_i] - snoice[t,:,:] = id.variables['snoice'][0,min_j:max_j,min_i:max_i] - meltt[t,:,:] = -1*id.variables['meltt'][0,min_j:max_j,min_i:max_i] - meltb[t,:,:] = -1*id.variables['meltb'][0,min_j:max_j,min_i:max_i] - meltl[t,:,:] = -1*id.variables['meltl'][0,min_j:max_j,min_i:max_i] - id.close() - t += 1 - - congel_avg = [] - frazil_avg = [] - snoice_avg = [] - meltt_avg = [] - meltb_avg = [] - meltl_avg = [] - for t in range(size(time)): - month = t % 12 - congel_avg.append(sum(congel[t,:,:]*aice[t,:,:]*dA)/sum(aice[t,:,:]*dA)*days_per_month[month]) - frazil_avg.append(sum(frazil[t,:,:]*aice[t,:,:]*dA)/sum(aice[t,:,:]*dA)*days_per_month[month]) - snoice_avg.append(sum(snoice[t,:,:]*aice[t,:,:]*dA)/sum(aice[t,:,:]*dA)*days_per_month[month]) - meltt_avg.append(sum(meltt[t,:,:]*aice[t,:,:]*dA)/sum(aice[t,:,:]*dA)*days_per_month[month]) - meltb_avg.append(sum(meltb[t,:,:]*aice[t,:,:]*dA)/sum(aice[t,:,:]*dA)*days_per_month[month]) - meltl_avg.append(sum(meltl[t,:,:]*aice[t,:,:]*dA)/sum(aice[t,:,:]*dA)*days_per_month[month]) - - congel_cum = cumsum(array(congel_avg)) - frazil_cum = cumsum(array(frazil_avg)) - snoice_cum = cumsum(array(snoice_avg)) - meltt_cum = cumsum(array(meltt_avg)) - meltb_cum = cumsum(array(meltb_avg)) - meltl_cum = cumsum(array(meltl_avg)) - total_cum = congel_cum + frazil_cum + snoice_cum + meltt_cum + meltb_cum + meltl_cum - - fig, ax = subplots(figsize=(8,6)) - ax.plot(time, congel_cum, label='Congelation', color='blue', linewidth=2) - ax.plot(time, frazil_cum, label='Frazil', color='red', linewidth=2) - ax.plot(time, snoice_cum, label='Snow-to-ice', color='cyan', linewidth=2) - ax.plot(time, meltt_cum, label='Top melt', color='magenta', linewidth=2) - ax.plot(time, meltb_cum, label='Basal melt', color='green', linewidth=2) - ax.plot(time, meltl_cum, label='Lateral melt', color='yellow', linewidth=2) - ax.plot(time, total_cum, label='Total', color='black', linewidth=2) - title_string = 'Cumulative volume tendency averaged over ' - if loc_flag == 'r': - title_string += 'Ross Sea' - elif loc_flag == 'w': - title_string += 'Weddell Sea' - title(title_string) - xlabel('Years') - ylabel('cm') - grid(True) - ax.legend(loc='upper left', prop=fontP) - if save: - fig.savefig(fig_name) - else: - fig.show() - - -if __name__ == "__main__": - - loc_flag = raw_input("Ross Sea (r) or Weddell Sea (w)? ") - action = raw_input("Save figure (s) or display on screen (d)? ") - if action == 's': - save = True - fig_name = raw_input("File name for figure: ") - else: - save = False - fig_name = None - nic_budget_ross_weddell(loc_flag, save, fig_name) - - - - diff --git a/paul_holland_hack.py b/paul_holland_hack.py index ee12947..3fbe700 100644 --- a/paul_holland_hack.py +++ b/paul_holland_hack.py @@ -2,25 +2,37 @@ from numpy import * from cartesian_grid_2d import * +# Paul Holland prevents spurious deep convection by applying an extra 1500 Gt/y +# of freshwater flux spread evenly south of 60S. Using the ROMS grid, calculate +# what the rate of input should be in kg/m^2/s. Print the value to the screen. +# Input: grid_file = path to ROMS grid file def paul_holland_hack (grid_file): total_fw = 1500 # Gt/y nbdry = -60 # Apply the freshwater evenly south of here sec_per_year = 365.25*24*60*60 + # Read grid and masks, making sure to get rid of the overlapping periodic + # boundary cells that are double-counted id = Dataset(grid_file, 'r') - lat = id.variables['lat_rho'][:,:] - lon = id.variables['lon_rho'][:,:] - mask_zice = id.variables['mask_zice'][:,:] - mask_rho = id.variables['mask_rho'][:,:] + lat = id.variables['lat_rho'][:,1:-1] + lon = id.variables['lon_rho'][:,1:-1] + mask_zice = id.variables['mask_zice'][:,1:-1] + mask_rho = id.variables['mask_rho'][:,1:-1] id.close() + # Mask out land and ice shelves mask = mask_rho - mask_zice + # Get differentials dx, dy = cartesian_grid_2d(lon, lat) + # Open ocean cells ocn_flag = mask == 1 + # Cells south of 60S loc_flag = lat < nbdry + # Total area of all open ocean cells total_area = sum(dx*dy*ocn_flag) print 'Total area = ' + str(total_area) + ' m^2' + # Total area of open ocean cells south of 60S target_area = sum(dx*dy*ocn_flag*loc_flag) print 'Area south of 60S = ' + str(target_area) + ' m^2' # Multiply by 1e12 to convert from Gt/y to kg/y @@ -30,6 +42,7 @@ def paul_holland_hack (grid_file): print 'Freshwater flux to add = ' + str(fw_flux) + 'kg/m^2/s' +# Command-line interface if __name__ == "__main__": grid_file = raw_input("Path to ROMS grid file: ") diff --git a/remove_iceshelves.py b/remove_iceshelves.py deleted file mode 100644 index b4e0de4..0000000 --- a/remove_iceshelves.py +++ /dev/null @@ -1,32 +0,0 @@ -from netCDF4 import Dataset -from numpy import * - -def remove_iceshelves (grid_file): - - id = Dataset(grid_file, 'a') - h = id.variables['h'][:,:] - mask_rho = id.variables['mask_rho'][:,:] - mask_zice = id.variables['mask_zice'][:,:] - - index = nonzero(mask_zice==1) - h[index] = 50 - mask_rho[index] = 0 - mask_u = mask_rho[:,1:]*mask_rho[:,:-1] - mask_v = mask_rho[1:,:]*mask_rho[:-1,:] - mask_psi = mask_rho[1:,1:]*mask_rho[1:,:-1]*mask_rho[:-1,1:]*mask_rho[:-1,:-1] - - id.variables['h'][:,:] = h - id.variables['mask_rho'][:,:] = mask_rho - id.variables['mask_u'][:,:] = mask_u - id.variables['mask_v'][:,:] = mask_v - id.variables['mask_psi'][:,:] = mask_psi - - id.close() - - -if __name__ == "__main__": - - grid_file = raw_input("Path to ROMS grid file to edit: ") - remove_iceshelves(grid_file) - - diff --git a/repeat_forcing.py b/repeat_forcing.py index 004c0a3..3b4b738 100644 --- a/repeat_forcing.py +++ b/repeat_forcing.py @@ -77,7 +77,7 @@ def process (directory, head, tail, year_start, perday, var_list): # User parameters to edit here # Data where files yyyy and yyyy exist - directory = '/short/m68/kaa561/ROMS-CICE-MCT/data/ERA_Interim/subdaily/30day_smoothed/' + directory = '/short/m68/kaa561/metroms_iceshelf/data/ERA_Interim/subdaily/30day_smoothed/' # First part of filename for AN and FC files an_head = 'AN_' fc_head = 'FC_' diff --git a/repeat_forcing_monthly.py b/repeat_forcing_monthly.py index a3fdc42..47d271d 100644 --- a/repeat_forcing_monthly.py +++ b/repeat_forcing_monthly.py @@ -61,7 +61,7 @@ def process (directory, head, tail, year_start): # User parameters to edit here # Data where files yyyy and yyyy exist - directory = '/short/m68/kaa561/ROMS-CICE-MCT/data/ERA_Interim/' + directory = '/short/m68/kaa561/metroms_iceshelf/data/ERA_Interim/' # First part of filename for AN and FC files an_head = 'AN_' fc_head = 'FC_' diff --git a/roms_ini_to_cice.py b/roms_ini_to_cice.py deleted file mode 100644 index fd71cc4..0000000 --- a/roms_ini_to_cice.py +++ /dev/null @@ -1,82 +0,0 @@ -import netCDF4 -import numpy as np - -def set_to_val(nc,variable,icemask,val=0.0,n=-999): - var = nc.variables[variable] - if n < 0: - var[:] = val*icemask - else: - var[n,:] = val*icemask - -# Make initialization of new cice-file... - -in_file = '/short/m68/kaa561/metroms_iceshelf/data/cice_ini_orig.nc' -cicefile = '/short/m68/kaa561/metroms_iceshelf/data/cice_ini.nc' - -NCAT = [0, 0.6445072, 1.391433, 2.470179, 4.567288, 1e+08] #upper limits of ice categories -varlist2d_null = ['uvel','vvel','scale_factor','swvdr','swvdf','swidr','swidf','strocnxT','strocnyT', - 'stressp_1','stressp_2','stressp_3','stressp_4','stressm_1','stressm_2','stressm_3', - 'stressm_4','stress12_1','stress12_2','stress12_3','stress12_4','iceumask','fsnow', 'accum_aice', 'accum_fresh', 'accum_fsalt', 'accum_fhocn', 'accum_fswthru', 'accum_strocnx', 'accum_strocny'] -varlist3d_null = ['vsnon','iage','alvl','vlvl','apnd','hpnd','ipnd','dhs','ffrac','qsno001'] - -nc_in = netCDF4.Dataset(in_file) -nc_cice = netCDF4.Dataset(cicefile,'r+') - -aice_r = nc_in.variables['aice'] # mean ice concentration for grid cell -aicen_c = nc_cice.variables['aicen'] # ice concentration per category in grid cell -vicen_c = nc_cice.variables['vicen'] # volume per unit area of ice (in category n) -tsfcn_c = nc_cice.variables['Tsfcn'] - -mask_r = nc_in.variables['mask'] - -hice_r = nc_in.variables['hice'] - -icemask = np.where((aice_r[:]*mask_r[:]) > 0.1, 1, 0) -print icemask - -for n in range(len(NCAT[1:])): - print 'Upper limit: '+str(NCAT[n+1]) - aicen_c[n,:,:] = np.where(np.logical_and( hice_r[:,:] > NCAT[n], hice_r[:,:] < NCAT[n+1] ), - aice_r[:,:],0.0) * mask_r[:,:] - vicen_c[n,:,:] = np.where(np.logical_and( hice_r[:,:] > NCAT[n], hice_r[:,:] < NCAT[n+1] ), - hice_r[:,:]*aice_r[:,:],0.0) * mask_r[:,:] - # Make sure no negative values exist: - aicen_c[n,:,:] = np.where(aicen_c[n,:,:] < 0.0, 0, aicen_c[n,:,:]) - vicen_c[n,:,:] = np.where(vicen_c[n,:,:] < 0.0, 0, vicen_c[n,:,:]) - tsfcn_c[n,:,:] = -20.15*icemask - 1.8*(1-icemask) - - for s in varlist3d_null: - print s - set_to_val(nc_cice,s,icemask,0.0,n) - set_to_val(nc_cice,'sice001',icemask,0.0,n) - set_to_val(nc_cice,'sice002',icemask,0.0,n) - set_to_val(nc_cice,'sice003',icemask,0.0,n) - set_to_val(nc_cice,'sice004',icemask,0.0,n) - set_to_val(nc_cice,'sice005',icemask,0.0,n) - set_to_val(nc_cice,'sice006',icemask,0.0,n) - set_to_val(nc_cice,'sice007',icemask,0.0,n) - set_to_val(nc_cice,'qice001',icemask,-2.0e8,n) - set_to_val(nc_cice,'qice002',icemask,-2.0e8,n) - set_to_val(nc_cice,'qice003',icemask,-2.0e8,n) - set_to_val(nc_cice,'qice004',icemask,-2.0e8,n) - set_to_val(nc_cice,'qice005',icemask,-2.0e8,n) - set_to_val(nc_cice,'qice006',icemask,-2.0e8,n) - set_to_val(nc_cice,'qice007',icemask,-2.0e8,n) - set_to_val(nc_cice,'alvl',icemask,1.0,n) - set_to_val(nc_cice,'vlvl',icemask,1.0,n) - -for s in varlist2d_null: - print s - set_to_val(nc_cice,s,icemask,0.0) - -nc_cice.istep1 = 0 -nc_cice.time = 0 -nc_cice.time_forc = 0 -nc_cice.nyr = 1 -nc_cice.month = 1 -nc_cice.mday = 1 -nc_cice.sec = 0 - - -nc_cice.close() -print 'done' diff --git a/romscice_atm_monthly.py b/romscice_atm_monthly.py index f2823a3..3e17b7f 100644 --- a/romscice_atm_monthly.py +++ b/romscice_atm_monthly.py @@ -21,11 +21,11 @@ def convert_file (year): # Paths of ROMS grid file, input ERA-Interim files, and output ROMS-CICE # files; other users will need to change these - grid_file = '/short/m68/kaa561/ROMS-CICE-MCT/apps/common/grid/circ30S_quarterdegree_10m.nc' + grid_file = '/short/m68/kaa561/metroms_iceshelf/apps/common/grid/circ30S_quarterdegree_10m.nc' input_atm_file = '/short/y99/kaa561/FESOM/ERA_Interim_monthly/AN_' + str(year) + '_monthly_orig.nc' input_ppt_file = '/short/y99/kaa561/FESOM/ERA_Interim_monthly/FC_' + str(year) + '_monthly_orig.nc' - output_atm_file = '/short/m68/kaa561/ROMS-CICE-MCT/data/ERA_Interim/monthly/historical/AN_' + str(year) + '_monthly.nc' - output_ppt_file = '/short/m68/kaa561/ROMS-CICE-MCT/data/ERA_Interim/monthly/historical/FC_' + str(year) + '_monthly.nc' + output_atm_file = '/short/m68/kaa561/metroms_iceshelf/data/ERA_Interim/monthly/historical/AN_' + str(year) + '_monthly.nc' + output_ppt_file = '/short/m68/kaa561/metroms_iceshelf/data/ERA_Interim/monthly/historical/FC_' + str(year) + '_monthly.nc' Lv = 2.5e6 # Latent heat of vapourisation, J/kg Rv = 461.5 # Ideal gas constant for water vapour, J/K/kg diff --git a/romscice_atm_subdaily.py b/romscice_atm_subdaily.py index 43f7943..ff51f43 100644 --- a/romscice_atm_subdaily.py +++ b/romscice_atm_subdaily.py @@ -32,11 +32,11 @@ def convert_file (year, count): # Paths of ROMS grid file, input ERA-Interim files, and output ROMS-CICE # files; other users will need to change these - grid_file = '/short/m68/kaa561/ROMS-CICE-MCT/apps/common/grid/circ30S_quarterdegree_good.nc' - input_atm_file = '/short/m68/kaa561/ROMS-CICE-MCT/data/ERA_Interim/AN_' + str(year) + '_subdaily_orig.nc' - input_ppt_file = '/short/m68/kaa561/ROMS-CICE-MCT/data/ERA_Interim/FC_' + str(year) + '_subdaily_orig.nc' - output_atm_file = '/short/m68/kaa561/ROMS-CICE-MCT/data/ERA_Interim/AN_' + str(year) + '_subdaily.nc' - output_ppt_file = '/short/m68/kaa561/ROMS-CICE-MCT/data/ERA_Interim/FC_' + str(year) + '_subdaily.nc' + grid_file = '/short/m68/kaa561/metroms_iceshelf/apps/common/grid/circ30S_quarterdegree_good.nc' + input_atm_file = '/short/m68/kaa561/metroms_iceshelf/data/ERA_Interim/AN_' + str(year) + '_subdaily_orig.nc' + input_ppt_file = '/short/m68/kaa561/metroms_iceshelf/data/ERA_Interim/FC_' + str(year) + '_subdaily_orig.nc' + output_atm_file = '/short/m68/kaa561/metroms_iceshelf/data/ERA_Interim/AN_' + str(year) + '_subdaily.nc' + output_ppt_file = '/short/m68/kaa561/metroms_iceshelf/data/ERA_Interim/FC_' + str(year) + '_subdaily.nc' logfile = str(year) + '.log' Lv = 2.5e6 # Latent heat of vapourisation, J/kg diff --git a/romscice_evap.py b/romscice_evap.py index e4346d4..e124e04 100644 --- a/romscice_evap.py +++ b/romscice_evap.py @@ -2,6 +2,16 @@ from numpy import * from interp_era2roms import * +# Read one year of ERA-Interim evaporation data (12-hourly), interpolate to the +# ROMS grid, and add to the existing FC forcing file created using +# romscice_atm_subdaily.nc. +# Input: year = integer containing the year to process +# count = time record in the given year to start with + +# This script only processes 50 12-hour timesteps at once to prevent memory +# overflow, and is designed to be called by a self-submitting batch script. +# See era_evap.job for an example. + def convert_file (year, count): # Make sure input arguments are integers (sometimes the batch script likes @@ -37,12 +47,13 @@ def convert_file (year, count): # Convert time units evap_time = evap_time/24.0 # days since 1900-01-01 00:00:0.0 evap_time = evap_time - 92*365 - 22 # days since 1992-01-01 00:00:0.0; note that there were 22 leap years between 1900 and 1992 - evap_time = evap_time - 0.5 # switch from precipitation over the preceding 12 hours to precipitation over the following 12 hours; this is easier for ROMS + evap_time = evap_time - 0.5 # switch from evaporation over the preceding 12 hours to evaporation over the following 12 hours; this is easier for ROMS # Also read ERA-Interim latitude and longitude lon_era = i_fid.variables['longitude'][:] lat_era = i_fid.variables['latitude'][:] i_fid.close() + # Define the variable in the output NetCDF file on the first timestep if count == 0: log = open(logfile, 'a') log.write('Adding evaporation to ' + output_evap_file + '\n') diff --git a/romscice_evap_era_monthly.py b/romscice_evap_era_monthly.py deleted file mode 100644 index 7aab620..0000000 --- a/romscice_evap_era_monthly.py +++ /dev/null @@ -1,231 +0,0 @@ -from numpy import * -from netCDF4 import Dataset, num2date -from matplotlib.pyplot import * - -def romscice_evap_era_monthly (roms_file, cice_file, month, colour_bounds=None, save=False, fig_name=None): - - era_head = '/short/y99/kaa561/CMIP5_forcing/atmos/climatology/ERA_Interim_monthly/ER_' - era_tail = '_monthly_orig.nc' - # Starting and ending days in each month - # Assume no leap years, we'll fix this later if we need - start_day = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1] - end_day = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31] - # Name of each month, for the title - month_name = ['January', 'February', 'March', 'April', 'May', 'June', 'July', 'August', 'September', 'October', 'November', 'December'] - # Degrees to radians conversion - deg2rad = pi/180.0 - # Conversion factor: m/12h to kg/m^2/s, opposite signs - era_conv = -1000./(12.*60.*60.) - # Conversion factor: cm/day to kg/m^2/s, opposite signs - cice_conv = -1000/(100.*24.*60.*60.) - - id_roms = Dataset(roms_file, 'r') - roms_lon = id_roms.variables['lon_rho'][1:-16,1:-1] - roms_lat = id_roms.variables['lat_rho'][1:-16,1:-1] - time_id = id_roms.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 - time = num2date(time_id[:], units=time_id.units, calendar=time_id.calendar.lower()) - - # Open the CICE file - id_cice = Dataset(cice_file, 'r') - - next_month = mod(month+1,12) - prev_month = mod(month-1,12) - - end_t = -1 - for t in range(size(time)-1, -1, -1): - if time[t].month-1 == month and time[t].day in range(end_day[month]-2, end_day[month]+1): - end_t = t - break - if time[t].month-1 == next_month and time[t].day in range(start_day[next_month], start_day[next_month]+2): - end_t = t - break - if end_t == -1: - print 'Error: ' + roms_file + ' does not contain a complete ' + month_name[month] - return - - start_t = -1 - for t in range(end_t, -1, -1): - if time[t].month-1 == prev_month and time[t].day in range(end_day[prev_month]-1, end_day[prev_month]+1): - start_t = t - break - if time[t].month-1 == month and time[t].day in range(start_day[month], start_day[month]+3): - start_t = t - break - if start_t == -1: - print 'Error: ' + roms_file + ' does not contain a complete ' + month_name[month] - return - - leap_year = False - roms_year = time[end_t].year - if month == 11: - roms_year = time[start_t].year - if mod(roms_year, 4) == 0: - leap_year = True - if mod(roms_year, 100) == 0: - leap_year = False - if mod(roms_year, 400) == 0: - leap_year = True - if leap_year: - end_day[1] = 29 - - romscice_evap = ma.empty(shape(roms_lon)) - romscice_evap[:,:] = 0.0 - num_days = 0 - - if time[start_t].month-1 == month and time[start_t].day == start_day[month]+2: - start_days = 5 - elif time[start_t].month-1 == month and time[start_t].day == start_day[month]+1: - start_days = 4 - elif time[start_t].month-1 == month and time[start_t].day == start_day[month]: - start_days = 3 - elif time[start_t].month-1 == prev_month and time[start_t].day == end_day[prev_month]: - start_days = 2 - elif time[start_t].month-1 == prev_month and time[start_t].day == end_day[prev_month]-1: - start_days = 1 - else: - print 'Error: starting index is month ' + str(time[start_t].month) + ', day ' + str(time[start_t].day) - return - - roms_evap = id_roms.variables['evaporation'][start_t,1:-16,1:-1] - cice_evap = id_cice.variables['evap_ai'][start_t,:-15,:]*cice_conv - aice = id_cice.variables['aice'][start_t,:-15,:] - romscice_evap += (aice*cice_evap + (1-aice)*roms_evap)*start_days - num_days += start_days - - for t in range(start_t+1, end_t): - roms_evap = id_roms.variables['evaporation'][t,1:-16,1:-1] - cice_evap = id_cice.variables['evap_ai'][t,:-15,:]*cice_conv - aice = id_cice.variables['aice'][t,:-15,:] - romscice_evap += (aice*cice_evap + (1-aice)*roms_evap)*5 - num_days += 5 - - if time[end_t].month-1 == next_month and time[end_t].day == start_day[next_month]+1: - end_days = 1 - elif time[end_t].month-1 == next_month and time[end_t].day == start_day[next_month]: - end_days = 2 - elif time[end_t].month-1 == month and time[end_t].day == end_day[month]: - end_days = 3 - elif time[end_t].month-1 == month and time[end_t].day == end_day[month]-1: - end_days = 4 - elif time[end_t].month-1 == month and time[end_t].day == end_day[month]-2: - end_days = 5 - else: - print 'Error: ending index is month ' + str(time[end_t].month) + ', day ' + str(time[end_t].day) - return - - roms_evap = id_roms.variables['evaporation'][end_t,1:-16,1:-1] - cice_evap = id_cice.variables['evap_ai'][end_t,:-15,:]*cice_conv - aice = id_cice.variables['aice'][end_t,:-15,:] - romscice_evap += (aice*cice_evap + (1-aice)*roms_evap)*end_days - num_days += end_days - - if num_days != end_day[month]: - print 'Error: found ' + str(num_days) + ' days instead of ' + str(end_day[month]) - return - - id_roms.close() - id_cice.close() - romscice_evap /= num_days - romscice_evap *= 1e6 - - id = Dataset(era_head + str(roms_year) + era_tail, 'r') - era_lon_1d = id.variables['longitude'][:] - era_lat_1d = id.variables['latitude'][:] - era_lon, era_lat = meshgrid(era_lon_1d, era_lat_1d) - era_evap = id.variables['e'][month,:,:]*era_conv - id.close() - era_evap *= 1e6 - - roms_x = -(roms_lat+90)*cos(roms_lon*deg2rad+pi/2) - roms_y = (roms_lat+90)*sin(roms_lon*deg2rad+pi/2) - era_x = -(era_lat+90)*cos(era_lon*deg2rad+pi/2) - era_y = (era_lat+90)*sin(era_lon*deg2rad+pi/2) - - if colour_bounds is not None: - lev = linspace(colour_bounds[0], colour_bounds[1], num=50) - else: - min_bound = min(amin(roms_evap), amin(era_evap)) - max_bound = max(amax(roms_evap), amax(era_evap)) - lev = linspace(min_bound, max_bound, num=50) - bdry = -30+90 - - fig = figure(figsize=(20,9)) - fig.add_subplot(1,2,1, aspect='equal') - contourf(era_x, era_y, era_evap, lev, extend='both') - title('ERA-Interim', fontsize=24) - xlim([-bdry, bdry]) - ylim([-bdry, bdry]) - axis('off') - fig.add_subplot(1,2,2, aspect='equal') - img = contourf(roms_x, roms_y, romscice_evap, lev, extend='both') - title('MetROMS', fontsize=24) - xlim([-bdry, bdry]) - ylim([-bdry, bdry]) - axis('off') - cbaxes = fig.add_axes([0.3, 0.04, 0.4, 0.04]) - cbar = colorbar(img, orientation='horizontal', cax=cbaxes) - suptitle(r'Evaporation (10$^{-6}$ kg/m$^2$/s), ' + month_name[month] + ' ' + str(roms_year), fontsize=30) - - if save: - fig.savefig(fig_name) - else: - fig.show() - - -if __name__ == "__main__": - - roms_file = raw_input("Path to ROMS file: ") - cice_file = raw_input("Path to CICE file with the same timesteps as ROMS file: ") - month = int(raw_input("Month number (1-12): "))-1 - colour_bounds = None - get_bounds = raw_input("Set bounds on colour scale (y/n)? ") - if get_bounds == 'y': - lower_bound = float(raw_input("Lower bound: ")) - upper_bound = float(raw_input("Upper bound: ")) - colour_bounds = [lower_bound, upper_bound] - action = raw_input("Save figure (s) or display on screen (d)? ") - if action == 's': - save = True - fig_name = raw_input("File name for figure: ") - elif action == 'd': - save = False - fig_name = None - romscice_evap_era_monthly(roms_file, cice_file, month, colour_bounds, save, fig_name) - - while True: - # Repeat until the user wants to exit - repeat = raw_input("Make another plot (y/n)? ") - if repeat == 'y': - while True: - # Ask for changes to the parameters until the user is done - changes = raw_input("Enter a parameter to change: (1) file paths, (2) month number, (3) colour bounds, (4) save/display; or enter to continue: ") - if len(changes) == 0: - # No more changes to parameters - break - else: - if int(changes) == 1: - # New ROMS file - roms_file = raw_input("Path to ROMS file: ") - cice_file = raw_input("Path to CICE file with the same timesteps as ROMS file: ") - elif int(changes) == 2: - # New month - month = int(raw_input("Month number (1-12): "))-1 - elif int(changes) == 3: - colour_bounds = None - get_bounds = raw_input("Set bounds on colour scale (y/n)? ") - if get_bounds == 'y': - lower_bound = float(raw_input("Lower bound: ")) - upper_bound = float(raw_input("Upper bound: ")) - colour_bounds = [lower_bound, upper_bound] - elif int(changes) == 4: - # Switch from save to display, or vice versa - save = not save - if save: - # Get a new figure name - fig_name = raw_input("File name for figure: ") - # Make the plot - romscice_evap_era_monthly(roms_file, cice_file, month, colour_bounds, save, fig_name) - else: - break diff --git a/romscice_gpcp.py b/romscice_gpcp.py index c5c2590..7af5123 100644 --- a/romscice_gpcp.py +++ b/romscice_gpcp.py @@ -2,18 +2,26 @@ from numpy import * from scipy.interpolate import RectBivariateSpline +# Convert one year of GPCP monthly averages for total precipitation to a +# ROMS-CICE input forcing file. +# Input: year = integer containing the year to process def romscice_gpcp (year): + # ROMS grid grid_file = '/short/m68/kaa561/metroms_iceshelf/apps/common/grid/circ30S_quarterdegree.nc' + # Beginning of GPCP file paths gpcp_head = '/short/m68/kaa561/gpcp/gpcp_cdr_v23rB1_y' + str(year) + '_m' + # Path to output file output_file = '/short/m68/kaa561/metroms_iceshelf/data/GPCP/precip_' + str(year) + '_monthly.nc' - conv_factor = 5e-4 # mm/day to m/12h + # Conversion factor: mm/day to m/12h + conv_factor = 5e-4 days_per_month = array([31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]) + # Check for leap years if year % 4 == 0: days_per_month[1] = 29 - # Read ROMS latitude and longitude + # Read ROMS grid id = Dataset(grid_file, 'r') lon_roms = id.variables['lon_rho'][:,:] lat_roms = id.variables['lat_rho'][:,:] @@ -83,6 +91,15 @@ def romscice_gpcp (year): out_id.close() +# Interpolate the given data from the GPCP grid to the ROMS grid. +# Input: +# A = nxm array containing data on the GPCP grid +# lon_gpcp = mx1 array containing GPCP longitude (0 to 360) +# lat_gpcp = nx1 array containing GPCP latitude +# lon_roms = pxq array containing ROMS longitude (0 to 360) +# lat_roms = pxq array containing ROMS latitude +# Output: +# B = pxq array containing data interpolated to the ROMS grid (lat x lon) def interp_gpcp2roms (A, lon_gpcp, lat_gpcp, lon_roms, lat_roms): # First build a function to approximate A with 2D splines diff --git a/romscice_ini_woa.py b/romscice_ini_woa.py index a3dfd8e..e53c9f7 100644 --- a/romscice_ini_woa.py +++ b/romscice_ini_woa.py @@ -264,12 +264,12 @@ def interp_woa2roms (A, lon_woa, lat_woa, depth_woa, lon_roms_3d, lat_roms_3d, z if __name__ == "__main__": # Path to ROMS grid file - grid_file = '../ROMS-CICE-MCT/apps/common/grid/circ30S_quarterdegree_good.nc' + grid_file = '../metroms_iceshelf/apps/common/grid/circ30S_quarterdegree.nc' # Path to World Ocean Atlas NetCDF file (converted from FESOM input using # woa_netcdf.py) woa_file = '/short/y99/kaa561/FESOM/woa01_ts.nc' # Path to desired output file - output_file = '../ROMS-CICE-MCT/data/woa_ini.nc' + output_file = '../metroms_iceshelf/data/woa_ini.nc' # Grid parameters: check grid_file and *.in file to make sure these are correct Tcline = 40 diff --git a/romscice_nbc.py b/romscice_nbc.py index d7da630..3307165 100644 --- a/romscice_nbc.py +++ b/romscice_nbc.py @@ -27,11 +27,11 @@ def convert_file (year): # Paths of ROMS grid file, input ECCO2 files (without the tail yyyymm.nc), # and output ROMS-CICE boundary condition file; other users will need to # change these - grid_file = '../ROMS-CICE-MCT/apps/common/grid/circ30S_quarterdegree_good.nc' - theta_base = '../ROMS-CICE-MCT/data/ECCO2/raw/THETA.1440x720x50.' + str(year) - salt_base = '../ROMS-CICE-MCT/data/ECCO2/raw/SALT.1440x720x50.' + str(year) - vvel_base = '../ROMS-CICE-MCT/data/ECCO2/raw/VVEL.1440x720x50.' + str(year) - output_file = '../ROMS-CICE-MCT/data/ECCO2/ecco2_cube92_lbc_' + str(year) + '.nc' + grid_file = '../metroms_iceshelf/apps/common/grid/circ30S_quarterdegree_good.nc' + theta_base = '../metroms_iceshelf/data/ECCO2/raw/THETA.1440x720x50.' + str(year) + salt_base = '../metroms_iceshelf/data/ECCO2/raw/SALT.1440x720x50.' + str(year) + vvel_base = '../metroms_iceshelf/data/ECCO2/raw/VVEL.1440x720x50.' + str(year) + output_file = '../metroms_iceshelf/data/ECCO2/ecco2_cube92_lbc_' + str(year) + '.nc' # Grid parameters; check grid_file and *.in to make sure these are correct Tcline = 40 diff --git a/romscice_nbc_rep.py b/romscice_nbc_rep.py index a03447c..bb3ec13 100644 --- a/romscice_nbc_rep.py +++ b/romscice_nbc_rep.py @@ -19,5 +19,5 @@ def run (filename): # User parameters if __name__ == "__main__": - filename = '../ROMS-CICE-MCT/data/ecco2_cube92_lbc_1995_rep.nc' + filename = '../metroms_iceshelf/data/ecco2_cube92_lbc_1995_rep.nc' run(filename) diff --git a/romscice_nbc_zeta.py b/romscice_nbc_zeta.py index 2557c0b..d06b724 100644 --- a/romscice_nbc_zeta.py +++ b/romscice_nbc_zeta.py @@ -17,8 +17,8 @@ def romscice_nbc_zeta (nbc_file): # Paths of ROMS grid file and input AVISO file - grid_file = '../ROMS-CICE-MCT/apps/common/grid/circ30S_quarterdegree_good.nc' - aviso_file = '../ROMS-CICE-MCT/data/originals/aviso_climatology.nc' + grid_file = '../metroms_iceshelf/apps/common/grid/circ30S_quarterdegree.nc' + aviso_file = '../metroms_iceshelf/data/originals/aviso_climatology.nc' # Read ROMS grid at northern boundary id = Dataset(grid_file, 'r') diff --git a/romscice_tides.py b/romscice_tides.py index a455579..3d9ab70 100644 --- a/romscice_tides.py +++ b/romscice_tides.py @@ -16,11 +16,11 @@ def romscice_tides (): # Path to ROMS grid file - grid_file = '../ROMS-CICE-MCT/apps/common/grid/circ30S_quarterdegree_10m.nc' + grid_file = '../metroms_iceshelf/apps/common/grid/circ30S_quarterdegree.nc' # Path to TPXO file - tpxo_file = '../ROMS-CICE-MCT/data/h_tpxo7.2.nc' + tpxo_file = '../metroms_iceshelf/data/h_tpxo7.2.nc' # Desired path to output file - out_file = '../ROMS-CICE-MCT/data/tides_tpxo72.nc' #_1yr.nc' + out_file = '../metroms_iceshelf/data/tides_tpxo72.nc' #_1yr.nc' # Bounds on latitude indices to read from TPXO2 (as close together as # possible while containing all the latitudes ROMS needs, so that there # aren't too many land points to fill with nearest neighbours which is diff --git a/seaice_budget_ross_weddell.py b/seaice_budget_ross_weddell.py deleted file mode 100644 index 70ac6b1..0000000 --- a/seaice_budget_ross_weddell.py +++ /dev/null @@ -1,136 +0,0 @@ -from netCDF4 import Dataset -from numpy import * -from matplotlib.pyplot import * -from matplotlib.font_manager import FontProperties -from cartesian_grid_2d import * - -def seaice_budget_ross_weddell (cice_file, loc_flag, save=False, fig_names=None): - - if loc_flag == 'r': - min_i = 650 - max_i = 850 - min_j = 80 - max_j = 200 - elif loc_flag == 'w': - min_i = 1150 - max_i = 1300 - min_j = 50 - max_j = 180 - - # Read CICE grid - id = Dataset(cice_file, 'r') - lon = id.variables['TLON'][:,:] - lat = id.variables['TLAT'][:,:] - # Calculate elements of area - dx, dy = cartesian_grid_2d(lon, lat) - dA = dx*dy - dA = dA[min_j:max_j, min_i:max_i] - # Read time values - time = id.variables['time'][:]/365.25 - # Read all the fields we need for timeseries - aice = id.variables['aice'][:,min_j:max_j,min_i:max_i] - dvidtt = id.variables['dvidtt'][:,min_j:max_j, min_i:max_i] - dvidtd = id.variables['dvidtd'][:,min_j:max_j, min_i:max_i] - congel = id.variables['congel'][:,min_j:max_j, min_i:max_i] - frazil = id.variables['frazil'][:,min_j:max_j, min_i:max_i] - snoice = id.variables['snoice'][:,min_j:max_j, min_i:max_i] - meltt = -1*id.variables['meltt'][:,min_j:max_j, min_i:max_i] - meltb = -1*id.variables['meltb'][:,min_j:max_j, min_i:max_i] - meltl = -1*id.variables['meltl'][:,min_j:max_j, min_i:max_i] - id.close() - - fontP = FontProperties() - fontP.set_size('small') - - dvidtt_avg = [] - dvidtd_avg = [] - for t in range(size(time)): - dvidtt_avg.append(sum(dvidtt[t,:,:]*aice[t,:,:]*dA)/sum(aice[t,:,:]*dA)) - dvidtd_avg.append(sum(dvidtd[t,:,:]*aice[t,:,:]*dA)/sum(aice[t,:,:]*dA)) - dvidtt_cum = cumsum(array(dvidtt_avg))*5 - dvidtd_cum = cumsum(array(dvidtd_avg))*5 - dvi_cum = dvidtt_cum + dvidtd_cum - - fig1, ax1 = subplots(figsize=(8,6)) - ax1.plot(time, dvidtt_cum, label='Thermodynamics', color='blue', linewidth=2) - ax1.plot(time, dvidtd_cum, label='Dynamics', color='green', linewidth=2) - ax1.plot(time, dvi_cum, label='Total', color='black', linewidth=2) - title_string = 'Cumulative volume tendency averaged over ' - if loc_flag == 'r': - title_string += 'Ross Sea' - elif loc_flag == 'w': - title_string += 'Weddell Sea' - title(title_string) - xlabel('Years') - ylabel('cm') - grid(True) - ax1.legend(loc='upper left', prop=fontP) - if save: - fig1.savefig(fig_names[0]) - else: - fig1.show() - - congel_avg = [] - frazil_avg = [] - snoice_avg = [] - meltt_avg = [] - meltb_avg = [] - meltl_avg = [] - for t in range(size(time)): - congel_avg.append(sum(congel[t,:,:]*aice[t,:,:]*dA)/sum(aice[t,:,:]*dA)) - frazil_avg.append(sum(frazil[t,:,:]*aice[t,:,:]*dA)/sum(aice[t,:,:]*dA)) - snoice_avg.append(sum(snoice[t,:,:]*aice[t,:,:]*dA)/sum(aice[t,:,:]*dA)) - meltt_avg.append(sum(meltt[t,:,:]*aice[t,:,:]*dA)/sum(aice[t,:,:]*dA)) - meltb_avg.append(sum(meltb[t,:,:]*aice[t,:,:]*dA)/sum(aice[t,:,:]*dA)) - meltl_avg.append(sum(meltl[t,:,:]*aice[t,:,:]*dA)/sum(aice[t,:,:]*dA)) - - congel_cum = cumsum(array(congel_avg))*5 - frazil_cum = cumsum(array(frazil_avg))*5 - snoice_cum = cumsum(array(snoice_avg))*5 - meltt_cum = cumsum(array(meltt_avg))*5 - meltb_cum = cumsum(array(meltb_avg))*5 - meltl_cum = cumsum(array(meltl_avg))*5 - total_cum = congel_cum + frazil_cum + snoice_cum + meltt_cum + meltb_cum + meltl_cum - - fig2, ax2 = subplots(figsize=(8,6)) - ax2.plot(time, congel_cum, label='Congelation', color='blue', linewidth=2) - ax2.plot(time, frazil_cum, label='Frazil', color='red', linewidth=2) - ax2.plot(time, snoice_cum, label='Snow-to-ice', color='cyan', linewidth=2) - ax2.plot(time, meltt_cum, label='Top melt', color='magenta', linewidth=2) - ax2.plot(time, meltb_cum, label='Basal melt', color='green', linewidth=2) - ax2.plot(time, meltl_cum, label='Lateral melt', color='yellow', linewidth=2) - ax2.plot(time, total_cum, label='Total', color='black', linewidth=2) - title_string = 'Cumulative volume tendency averaged over ' - if loc_flag == 'r': - title_string += 'Ross Sea' - elif loc_flag == 'w': - title_string += 'Weddell Sea' - title(title_string) - xlabel('Years') - ylabel('cm') - grid(True) - ax2.legend(loc='lower left', prop=fontP) - if save: - fig2.savefig(fig_names[1]) - else: - fig2.show() - - -if __name__ == "__main__": - - cice_file = raw_input("Path to CICE history file: ") - loc_flag = raw_input("Ross Sea (r) or Weddell Sea (w)? ") - action = raw_input("Save figures (s) or display on screen (d)? ") - if action == 's': - save = True - name1 = raw_input("File name for first figure (thermodynamic versus dynamic): ") - name2 = raw_input("File name for second figure (thermodynamic terms): ") - fig_names = [name1, name2] - else: - save = False - fig_names = None - seaice_budget_ross_weddell(cice_file, loc_flag, save, fig_names) - - - - diff --git a/seaice_budget_thickice.py b/seaice_budget_thickice.py deleted file mode 100644 index e89343f..0000000 --- a/seaice_budget_thickice.py +++ /dev/null @@ -1,119 +0,0 @@ -from netCDF4 import Dataset -from numpy import * -from matplotlib.pyplot import * -from matplotlib.font_manager import FontProperties -from cartesian_grid_2d import * - -def seaice_budget_thickice (cice_file, save=False, fig_names=None): - - threshold = 5.0 - - # Read CICE grid - id = Dataset(cice_file, 'r') - lon = id.variables['TLON'][:,:] - lat = id.variables['TLAT'][:,:] - # Calculate elements of area - dx, dy = cartesian_grid_2d(lon, lat) - dA = dx*dy - # Read time values - time = id.variables['time'][:]/365.25 - # Read final sea ice thickness - hi_f = id.variables['hi'][-1,:,:] - # Read all the fields we need for timeseries - aice = id.variables['aice'][:,:,:] - dvidtt = id.variables['dvidtt'][:,:,:] - dvidtd = id.variables['dvidtd'][:,:,:] - congel = id.variables['congel'][:,:,:] - frazil = id.variables['frazil'][:,:,:] - snoice = id.variables['snoice'][:,:,:] - meltt = -1*id.variables['meltt'][:,:,:] - meltb = -1*id.variables['meltb'][:,:,:] - meltl = -1*id.variables['meltl'][:,:,:] - id.close() - - flag = hi_f > threshold - - fontP = FontProperties() - fontP.set_size('small') - - dvidtt_avg = [] - dvidtd_avg = [] - for t in range(size(time)): - dvidtt_avg.append(sum(dvidtt[t,:,:]*aice[t,:,:]*dA*flag)/sum(aice[t,:,:]*dA*flag)) - dvidtd_avg.append(sum(dvidtd[t,:,:]*aice[t,:,:]*dA*flag)/sum(aice[t,:,:]*dA*flag)) - dvidtt_cum = cumsum(array(dvidtt_avg))*5 - dvidtd_cum = cumsum(array(dvidtd_avg))*5 - dvi_cum = dvidtt_cum + dvidtd_cum - - fig1, ax1 = subplots(figsize=(8,6)) - ax1.plot(time, dvidtt_cum, label='Thermodynamics', color='blue', linewidth=2) - ax1.plot(time, dvidtd_cum, label='Dynamics', color='green', linewidth=2) - ax1.plot(time, dvi_cum, label='Total', color='black', linewidth=2) - title('Cumulative volume tendency averaged over thick ice regions') - xlabel('Years') - ylabel('cm') - grid(True) - ax1.legend(loc='upper left', prop=fontP) - if save: - fig1.savefig(fig_names[0]) - else: - fig1.show() - - congel_avg = [] - frazil_avg = [] - snoice_avg = [] - meltt_avg = [] - meltb_avg = [] - meltl_avg = [] - for t in range(size(time)): - congel_avg.append(sum(congel[t,:,:]*aice[t,:,:]*dA*flag)/sum(aice[t,:,:]*dA*flag)) - frazil_avg.append(sum(frazil[t,:,:]*aice[t,:,:]*dA*flag)/sum(aice[t,:,:]*dA*flag)) - snoice_avg.append(sum(snoice[t,:,:]*aice[t,:,:]*dA*flag)/sum(aice[t,:,:]*dA*flag)) - meltt_avg.append(sum(meltt[t,:,:]*aice[t,:,:]*dA*flag)/sum(aice[t,:,:]*dA*flag)) - meltb_avg.append(sum(meltb[t,:,:]*aice[t,:,:]*dA*flag)/sum(aice[t,:,:]*dA*flag)) - meltl_avg.append(sum(meltl[t,:,:]*aice[t,:,:]*dA*flag)/sum(aice[t,:,:]*dA*flag)) - - congel_cum = cumsum(array(congel_avg)) - frazil_cum = cumsum(array(frazil_avg)) - snoice_cum = cumsum(array(snoice_avg)) - meltt_cum = cumsum(array(meltt_avg)) - meltb_cum = cumsum(array(meltb_avg)) - meltl_cum = cumsum(array(meltl_avg)) - total_cum = congel_cum + frazil_cum + snoice_cum + meltt_cum + meltb_cum + meltl_cum - - fig2, ax2 = subplots(figsize=(8,6)) - ax2.plot(time, congel_cum, label='Congelation', color='blue', linewidth=2) - ax2.plot(time, frazil_cum, label='Frazil', color='red', linewidth=2) - ax2.plot(time, snoice_cum, label='Snow-to-ice', color='cyan', linewidth=2) - ax2.plot(time, meltt_cum, label='Top melt', color='magenta', linewidth=2) - ax2.plot(time, meltb_cum, label='Basal melt', color='green', linewidth=2) - ax2.plot(time, meltl_cum, label='Lateral melt', color='yellow', linewidth=2) - ax2.plot(time, total_cum, label='Total', color='black', linewidth=2) - title('Cumulative volume tendency averaged over thick ice regions') - xlabel('Years') - ylabel('cm') - grid(True) - ax2.legend(loc='lower left', prop=fontP) - if save: - fig2.savefig(fig_names[1]) - else: - fig2.show() - - -if __name__ == "__main__": - - cice_file = raw_input("Path to CICE history file: ") - action = raw_input("Save figures (s) or display on screen (d)? ") - if action == 's': - save = True - name1 = raw_input("File name for first figure (thermodynamic versus dynamic): ") - name2 = raw_input("File name for second figure (thermodynamic terms): ") - fig_names = [name1, name2] - else: - save = False - fig_names = None - seaice_budget_thickice(cice_file, save, fig_names) - - - - diff --git a/seasonal_avg_cice.py b/seasonal_avg_cice.py index 58b1af5..69830fb 100644 --- a/seasonal_avg_cice.py +++ b/seasonal_avg_cice.py @@ -1,6 +1,16 @@ from netCDF4 import Dataset, num2date from numpy import * +# Calculate seasonal averages (DJF, MAM, JJA, SON) of the given variable in the +# given CICE file. +# Input: +# file_path = path to CICE output file containing 5-day averages, including at +# least one complete December-November period. If there are multiple +# such instances the last one will be plotted. +# var = variable name +# shape = vector containing the dimensions (excluding time) of the variable +# Output: +# seasonal_data = array of data averaged over each season, dimension 4 x shape def seasonal_avg_cice (file_path, var, shape): # Starting and ending months (1-based) for each season diff --git a/seasonal_avg_roms.py b/seasonal_avg_roms.py index c3d5eca..efa959e 100644 --- a/seasonal_avg_roms.py +++ b/seasonal_avg_roms.py @@ -1,6 +1,16 @@ from netCDF4 import Dataset, num2date from numpy import * +# Calculate seasonal averages (DJF, MAM, JJA, SON) of the given variable in the +# given ROMS file. +# Input: +# file_path = path to ROMS output file containing 5-day averages, including at +# least one complete December-November period. If there are multiple +# such instances the last one will be plotted. +# var = variable name +# shape = vector containing the dimensions (excluding time) of the variable +# Output: +# seasonal_data = array of data averaged over each season, dimension 4 x shape def seasonal_avg_roms (file_path, var, shape): # Starting and ending months (1-based) for each season diff --git a/snow_balance.py b/snow_balance.py deleted file mode 100644 index 1cf3bc9..0000000 --- a/snow_balance.py +++ /dev/null @@ -1,53 +0,0 @@ -from netCDF4 import Dataset -from numpy import * -from matplotlib.pyplot import * - -def snow_balance (cice_file): - - deg2rad = pi/180.0 - - id = Dataset(cice_file, 'r') - lon_tmp = id.variables['TLON'][:-15,:] - lat_tmp = id.variables['TLAT'][:-15,:] - hs_i = id.variables['hs'][0,:-15,:] - hs_f = id.variables['hs'][72,:-15,:] - total_snow = sum(id.variables['snow_ai'][0:73,:-15,:]*5, axis=0)/100.0 - total_snoice = sum(id.variables['snoice'][0:73,:-15,:]*5, axis=0)/100.0 - total_melts = sum(id.variables['melts'][0:73,:-15,:]*5, axis=0)/100.0 - id.close() - - balance_tmp = hs_i + total_snow - total_snoice - total_melts - hs_f - - lon = ma.empty([size(lon_tmp,0), size(lon_tmp,1)+1]) - lat = ma.empty([size(lat_tmp,0), size(lat_tmp,1)+1]) - balance = ma.empty([size(balance_tmp,0), size(balance_tmp,1)+1]) - lon[:,:-1] = lon_tmp - lon[:,-1] = lon_tmp[:,0] - lat[:,:-1] = lat_tmp - lat[:,-1] = lat_tmp[:,0] - balance[:,:-1] = balance_tmp - balance[:,-1] = balance_tmp[:,0] - - x = -(lat+90)*cos(lon*deg2rad+pi/2) - y = (lat+90)*sin(lon*deg2rad+pi/2) - lev = linspace(-2, 0, num=50) - - fig = figure(figsize=(16,12)) - fig.add_subplot(1,1,1, aspect='equal') - contourf(x, y, balance, lev, extend='both') - cbar = colorbar() - cbar.ax.tick_params(labelsize=20) - title('Snow balance (m) over first year of simulation\nTotal snow_ai - total melts - total snoice - change in hs', fontsize=30) - axis('off') - - fig.show() - #fig.savefig('snow_balance.png') - - -if __name__ == "__main__": - - cice_file = raw_input("Path to CICE history file, containing at least one year of 5-day averages: ") - snow_balance(cice_file) - - - diff --git a/snow_budget.py b/snow_budget.py deleted file mode 100644 index 6adcb97..0000000 --- a/snow_budget.py +++ /dev/null @@ -1,61 +0,0 @@ -from netCDF4 import Dataset -from numpy import * -from matplotlib.pyplot import * -from cartesian_grid_2d import * - -def snow_budget (cice_file): - - # Read CICE grid - id = Dataset(cice_file, 'r') - lon = id.variables['TLON'][:,:] - lat = id.variables['TLAT'][:,:] - # Calculate elements of area - dx, dy = cartesian_grid_2d(lon, lat) - dA = dx*dy - # Read time values - time = id.variables['time'][:]/365.25 - hs = id.variables['hs'][:,:,:] - aice = id.variables['aice'][:,:,:] - snow_ai = id.variables['snow_ai'][:,:,:]/100 - snoice = -1*id.variables['snoice'][:,:,:]/100/3.6 - melts = -1*id.variables['melts'][:,:,:]/100/3.6 - id.close() - - vs = [] - total_snow = [] - total_snoice = [] - total_melts = [] - for t in range(size(time)): - vs.append(sum(hs[t,:,:]*aice[t,:,:]*dA)) - total_snow.append(sum(snow_ai[t,:,:]*dA)) - total_snoice.append(sum(snoice[t,:,:]*aice[t,:,:]*dA)) - total_melts.append(sum(melts[t,:,:]*aice[t,:,:]*dA)) - delta_vs = vs - vs[0] - cum_snow = cumsum(array(total_snow))*5 - cum_snoice = cumsum(array(total_snoice))*5 - cum_melts = cumsum(array(total_melts))*5 - total = cum_snow + cum_snoice + cum_melts - - fig, ax = subplots(figsize=(12,9)) - ax.plot(time, delta_vs/1e9, label='Snow volume', color='black', linewidth=2) - ax.plot(time, cum_snow/1e9, label='Total snowfall', color='red', linewidth=2) - ax.plot(time, cum_snoice/1e9, label='Total snow-to-ice', color='blue', linewidth=2) - ax.plot(time, cum_melts/1e9, label='Total snow melt', color='yellow', linewidth=2) - ax.plot(time, total/1e9, label='Expected snow volume', color='green', linewidth=2) - xlabel('Time (years)') - ylabel(r'km$^3$') - grid(True) - ax.legend(loc='lower left') - fig.show() - - -if __name__ == "__main__": - - cice_file = raw_input("Path to CICE history file: ") - snow_budget(cice_file) - - - - - - diff --git a/ssflux_monthly.py b/ssflux_monthly.py deleted file mode 100644 index 23fec3d..0000000 --- a/ssflux_monthly.py +++ /dev/null @@ -1,91 +0,0 @@ -from netCDF4 import Dataset, num2date -from numpy import * -from matplotlib.pyplot import * -from monthly_avg_roms import * - -def ssflux_monthly (roms_file, month, bound=None, save=False, fig_name=None): - - # Month names for plot titles - month_name = ['January', 'February', 'March', 'April', 'May', 'June', 'July', 'August', 'September', 'October', 'November', 'December'] - deg2rad = pi/180 - - # Read ROMS grid - id = Dataset(roms_file, 'r') - lon = id.variables['lon_rho'][:-15,:-1] - lat = id.variables['lat_rho'][:-15,:-1] - num_lon = id.variables['lon_rho'].shape[1] - num_lat = id.variables['lat_rho'].shape[0] - id.close() - - ssflux = monthly_avg_roms(roms_file, 'ssflux', [num_lat, num_lon], month) - ssflux = ssflux[:-15,:-1] - ssflux *= 1e6 - - x = -(lat+90)*cos(lon*deg2rad+pi/2) - y = (lat+90)*sin(lon*deg2rad+pi/2) - - if bound is None: - bound = amax(abs(ssflux)) - lev = linspace(-bound, bound, num=50) - - fig = figure(figsize=(16,12)) - fig.add_subplot(1,1,1,aspect='equal') - img = contourf(x, y, ssflux, lev, cmap='RdYlBu_r', extend='both') - cbar = colorbar() - cbar.ax.tick_params(labelsize=20) - title(month_name[month] + r' surface salinity flux (10$^{-6}$ kg/m$^2$/s)', fontsize=30) - axis('off') - - if save: - fig.savefig(fig_name) - else: - fig.show() - - -if __name__ == "__main__": - - roms_file = raw_input("Path to ROMS output file: ") - month = int(raw_input("Month number (1-12): "))-1 - bound = None - get_bound = raw_input("Set bounds on colour scale (y/n)? ") - if get_bound == 'y': - bound = float(raw_input("Maximum absolute value (1e-6 kg/m^2/s): ")) - 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 - ssflux_monthly(roms_file, month, bound, save, fig_name) - - while True: - repeat = raw_input("Make another plot (y/n)? ") - if repeat == 'y': - while True: - changes = raw_input("Enter a parameter to change: (1) file path, (2) month, (3) colour bounds, (4) save/display; or enter to continue: ") - if len(changes) == 0: - break - else: - if int(changes) == 1: - roms_file = raw_input("Path to ROMS output file: ") - elif int(changes) == 2: - month = int(raw_input("Month number (1-12): "))-1 - elif int(changes) == 3: - get_bound = raw_input("Set bounds on colour scale (y/n)? ") - if get_bound == 'y': - bound = float(raw_input("Maximum absolute value (1e-6 kg/m^2/s): ")) - elif int(changes) == 4: - save = not save - if save: - fig_name = raw_input("File name for figure: ") - ssflux_monthly(roms_file, month, bound, save, fig_name) - else: - break - - - - - - - diff --git a/ssflux_tamura_annual.py b/ssflux_tamura_annual.py index 04c8eb3..209bf1d 100644 --- a/ssflux_tamura_annual.py +++ b/ssflux_tamura_annual.py @@ -2,8 +2,17 @@ from netCDF4 import Dataset, num2date from matplotlib.pyplot import * +# Make a figure comparing annually-averaged sea ice to ocean salt fluxes, +# from Tamura's dataset to CICE output. +# Input: +# cice_file = path to annually averaged CICE file +# year = year of Tamura data to plot +# save = optional boolean to save the figure to a file, rather than displaying +# it on the screen +# fig_name = if save=True, path to the desired filename for figure def ssflux_tamura_annual (cice_file, year, save=False, fig_name=None): + # Path to Tamura file for this year tamura_file = '/short/m68/kaa561/tamura_fluxes/Tamura_ssflux_' + str(year) + '_monthly.nc' # Degrees to radians conversion deg2rad = pi/180.0 diff --git a/temp_salt_slice.py b/temp_salt_slice.py index ae2808e..87f8655 100644 --- a/temp_salt_slice.py +++ b/temp_salt_slice.py @@ -4,6 +4,16 @@ from calc_z import * from interp_lon_roms import * +# Create a 2x1 plot showing zonal slices (depth vs latitude) of temperature and +# salinity interpolated to the given longitude, at the given timestep. +# Input: +# file_path = path to ROMS history or averages file +# tstep = time index in file_path to plot (1-indexed) +# lon0 = longitude to interpolate to (-180 to 180) +# depth_min = deepest depth to plot (negative, in metres) +# save = optional boolean flag; if True, the figure will be saved with file name +# fig_name, if False, the figure will display on the screen +# fig_name = optional string containing filename for figure, if save=True def temp_salt_slice (file_path, tstep, lon0, depth_min, save=False, fig_name=None): # Grid parameters @@ -12,12 +22,15 @@ def temp_salt_slice (file_path, tstep, lon0, depth_min, save=False, fig_name=Non hc = 40 N = 31 + # Month names for titles month_names = ['January', 'February', 'March', 'April', 'May', 'June', 'July', 'August', 'September', 'October', 'November', 'December'] + # Bounds on colour scales for temperature and salinity var_min = [-2, 33.8] var_max = [3, 34.8] var_tick = [1, 0.2] + # Read temperature, salinity, and grid variables id = Dataset(file_path, 'r') temp_3d = id.variables['temp'][tstep-1,:,:-15,:] salt_3d = id.variables['salt'][tstep-1,:,:-15,:] @@ -26,6 +39,7 @@ def temp_salt_slice (file_path, tstep, lon0, depth_min, save=False, fig_name=Non zice = id.variables['zice'][:-15,:] lon_2d = id.variables['lon_rho'][:-15,:] lat_2d = id.variables['lat_rho'][:-15,:] + # Read time axis and convert to Date objects time_id = id.variables['ocean_time'] time = num2date(time_id[tstep-1], units=time_id.units, calendar=time_id.calendar.lower()) id.close() @@ -33,16 +47,20 @@ def temp_salt_slice (file_path, tstep, lon0, depth_min, save=False, fig_name=Non # 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) + # Get the date for the title date_string = str(time.day) + ' ' + month_names[time.month-1] + ' ' + str(time.year) + # 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' + # Make sure we are in the range 0-360 if lon0 < 0: lon0 += 360 + # Interpolate temperature, salinity, z, and latitude to lon0 temp, z, lat = interp_lon_roms(temp_3d, z_3d, lat_2d, lon_2d, lon0) salt, z, lat = interp_lon_roms(salt_3d, z_3d, lat_2d, lon_2d, lon0) @@ -69,10 +87,13 @@ def temp_salt_slice (file_path, tstep, lon0, depth_min, save=False, fig_name=Non # Show the first 2 degrees of this land mask lat_max = max(lat[:,j_max]) + 2 + # Colour levels lev1 = linspace(var_min[0], var_max[0], num=50) lev2 = linspace(var_min[1], var_max[1], num=50) + # Plot fig = figure(figsize=(24,6)) + # Temperature ax = fig.add_subplot(1,2,1) img1 = contourf(lat, z, temp, lev1, extend='both') xlim([lat_min, lat_max]) @@ -82,6 +103,7 @@ def temp_salt_slice (file_path, tstep, lon0, depth_min, save=False, fig_name=Non title(r'Temperature ($^{\circ}$C)', fontsize=20) cbar1 = colorbar(img1, ticks=arange(var_min[0], var_max[0]+var_tick[0], var_tick[0])) cbar1.ax.tick_params(labelsize=16) + # Salinity ax = fig.add_subplot(1,2,2) img2 = contourf(lat, z, salt, lev2, extend='both') xlim([lat_min, lat_max]) @@ -94,15 +116,18 @@ def temp_salt_slice (file_path, tstep, lon0, depth_min, save=False, fig_name=Non suptitle(date_string + ', ' + lon_string, fontsize=24) subplots_adjust(wspace=0.025) + # Finished if save: fig.savefig(fig_name) else: fig.show() + # Convert back to the range -180 to 180, in case this script repeats if lon0 > 180: lon0 -= 360 +# Command-line interface if __name__ == "__main__": file_path = raw_input("Path to ocean averages file: ") @@ -116,14 +141,18 @@ def temp_salt_slice (file_path, tstep, lon0, depth_min, save=False, fig_name=Non elif action == 'd': save = False fig_name = None + # Make the plot temp_salt_slice(file_path, tstep, lon0, depth_min, save, fig_name) + # Repeat until the user wants to exit while True: repeat = raw_input("Make another plot (y/n)? ") if repeat == 'y': + # Ask for changes to parameters until the user is done while True: changes = raw_input("Enter a parameter to change: (1) file path, (2) time index, (3) longitude, (4) deepest depth, (5) save/display; or enter to continue: ") if len(changes) == 0: + # No more changes to parameters break else: if int(changes) == 1: @@ -137,9 +166,12 @@ def temp_salt_slice (file_path, tstep, lon0, depth_min, save=False, fig_name=Non elif int(changes) == 5: save = not save if save: + # Get new figure name fig_name = raw_input("File name for figure: ") + # Make the plot temp_salt_slice(file_path, tstep, lon0, depth_min, save, fig_name) else: + # No more figures break diff --git a/timeseries_i2osalt.py b/timeseries_i2osalt.py index c8f4e62..4bb3265 100644 --- a/timeseries_i2osalt.py +++ b/timeseries_i2osalt.py @@ -4,14 +4,25 @@ from os.path import * from cartesian_grid_2d import * +# Calculate and plot timeseries of area-averaged sea ice to ocean salt flux +# during a ROMS-CICE simulation. +# Input: +# file_path = path to CICE history file +# log_path = path to log file (if it exists, previously calculated values will +# be read from it; regardless, it will be overwritten with all +# calculated values following computation) def timeseries_i2osalt (file_path, log_path): + # Density of freshwater rho_fw = 1000.0 - rho_sw = 1026.0 + # Density of seawater + rho_sw = 1025.0 + # Conversion from m/s to cm/day mps_to_cmpday = 8.64e6 time = [] avg_ssflux = [] + # Check if the log file exists if exists(log_path): print 'Reading previously calculated values' f = open(log_path, 'r') @@ -41,13 +52,20 @@ def timeseries_i2osalt (file_path, log_path): time.append(new_time[t]) print 'Reading data' + # Read freshwater, salt, and rain fluxes (all scaled by aice) and + # sea surface salinity + # Throw away northern sponge layer fresh_ai = id.variables['fresh_ai'][:,:-15,:] sss = id.variables['sss'][:,:-15,:] rain_ai = id.variables['rain_ai'][:,:-15,:] fsalt_ai = id.variables['fsalt_ai'][:,:-15,:] id.close() + # Build timeseries for t in range(size(new_time)): + # Merge CICE's freshwater and salt fluxes as in set_vbc.F + # Subtract rain because we don't care about that + # Convert to kg/m^2/s avg_ssflux.append(sum(-1/rho_fw*((fresh_ai[t,:,:]-rain_ai[t,:,:])*sss[t,:,:]*rho_sw/mps_to_cmpday - fsalt_ai[t,:,:]*1e3)*dA)/sum(dA)) print 'Plotting' @@ -69,6 +87,7 @@ def timeseries_i2osalt (file_path, log_path): f.close() +# Command-line interface if __name__ == "__main__": file_path = raw_input("Path to CICE history file: ") diff --git a/timeseries_sss.py b/timeseries_sss.py index dc4e2b3..899b5eb 100644 --- a/timeseries_sss.py +++ b/timeseries_sss.py @@ -4,11 +4,19 @@ from os.path import * from cartesian_grid_2d import * +# Calculate and plot timeseries of area-averaged sea surface salinity and +# surface salt flux during a ROMS simulation. +# Input: +# file_path = path to ROMS averages file +# log_path = path to log file (if it exists, previously calculated values will +# be read from it; regardless, it will be overwritten with all +# calculated values following computation) def timeseries_sss (file_path, log_path): time = [] avg_sss = [] avg_ssflux = [] + # Check if the log file exists if exists(log_path): print 'Reading previously calculated values' f = open(log_path, 'r') @@ -34,10 +42,9 @@ def timeseries_sss (file_path, log_path): lon = id.variables['lon_rho'][:-15,1:-1] lat = id.variables['lat_rho'][:-15,1:-1] zice = id.variables['zice'][:-15,1:-1] - + # Calculate area on the tracer grid dx, dy = cartesian_grid_2d(lon, lat) dA = ma.masked_where(zice!=0, dx*dy) - # Read time values and convert from seconds to years new_time = id.variables['ocean_time'][:]/(365.25*24*60*60) # Concatenate with time values from log file @@ -45,10 +52,13 @@ def timeseries_sss (file_path, log_path): time.append(new_time[t]) print 'Reading data' + # Read surface salinity and salt flux + # Throw away overlapping periodic boundary and northern sponge layer sss = id.variables['salt'][:,-1,:-15,1:-1] ssflux = id.variables['ssflux'][:,:-15,1:-1] id.close() + # Build timeseries for t in range(size(new_time)): avg_sss.append(sum(sss[t,:,:]*dA)/sum(dA)) avg_ssflux.append(sum(ssflux[t,:,:]*dA)/sum(dA)) @@ -82,6 +92,7 @@ def timeseries_sss (file_path, log_path): f.close() +# Command-line interface if __name__ == "__main__": file_path = raw_input("Path to ocean history/averages file: ") diff --git a/woa_jan_netcdf.py b/woa_jan_netcdf.py deleted file mode 100644 index b7ad920..0000000 --- a/woa_jan_netcdf.py +++ /dev/null @@ -1,95 +0,0 @@ -from netCDF4 import Dataset -from numpy import * - -def woa_jan_netcdf (): - - temp_file = '/short/m68/kaa561/woa_jan/woa13_decav_t01mn01v2.csv' - salt_file = '/short/m68/kaa561/woa_jan/woa13_decav_t01mn01v2.csv' - out_file = '/short/m68/kaa561/woa_jan/woa_jan.nc' - - lat = arange(-77.5, 87.5+1, 1) - lon = arange(-179.5, 179.5+1, 1) - num_lat = size(lat) - num_lon = size(lon) - - file = open(temp_file, 'r') - file.readline() - depth_data = file.readline() - index = depth_data.index('0') - depth_data = depth_data[index:] - depth_vals = depth_data.split(',') - depth = [] - for elm in depth_vals: - depth.append(float(elm)) - depth = array(depth) - num_depth = size(depth) - - temp_data = zeros([num_depth, num_lat, num_lon]) - temp_data[:,:,:] = -999 - for line in file: - data = line.split(',') - lat0 = float(data[0]) - lon0 = float(data[1]) - j = where(lat==lat0)[0][0] - i = where(lon==lon0)[0][0] - data = data[2:] - k = 0 - for elm in data: - if elm not in ['', '\n']: - temp_data[k,j,i] = float(elm) - k += 1 - file.close() - temp_data = ma.masked_where(temp_data==-999, temp_data) - - salt_data = zeros([num_depth, num_lat, num_lon]) - salt_data[:,:,:] = -999 - file = open(salt_file, 'r') - file.readline() - file.readline() - for line in file: - data = line.split(',') - lat0 = float(data[0]) - lon0 = float(data[1]) - j = where(lat==lat0)[0][0] - i = where(lon==lon0)[0][0] - data = data[2:] - k = 0 - for elm in data: - if elm not in ['', '\n']: - salt_data[k,j,i] = float(elm) - k += 1 - file.close() - salt_data = ma.masked_where(salt_data==-999, salt_data) - - id = Dataset(out_file, 'w') - id.createDimension('longitude', num_lon) - id.createDimension('latitude', num_lat) - id.createDimension('depth', num_depth) - id.createVariable('longitude', 'f8', ('longitude')) - id.variables['longitude'].units = 'degrees' - id.variables['longitude'][:] = lon - id.createVariable('latitude', 'f8', ('latitude')) - id.variables['latitude'].units = 'degrees' - id.variables['latitude'][:] = lat - id.createVariable('depth', 'f8', ('depth')) - id.variables['depth'].units = 'metres' - id.variables['depth'][:] = depth - id.createVariable('temp', 'f8', ('depth', 'latitude', 'longitude')) - id.variables['temp'].units = 'C' - id.variables['temp'][:,:,:] = temp_data - id.createVariable('salt', 'f8', ('depth', 'latitude', 'longitude')) - id.variables['salt'].units = 'psu' - id.variables['salt'][:,:,:] = salt_data - id.close() - - -if __name__ == "__main__": - - woa_jan_netcdf() - - - - - - -