From 75bb31b28cb988ead90c60aea94a39f89c145999 Mon Sep 17 00:00:00 2001 From: Kaitlin Alexander Date: Tue, 11 Oct 2016 11:49:23 +1100 Subject: [PATCH] Moved CMIP5 code to separate repository; added i-slice plots --- .gitignore | 6 + aice_hi_seasonal.py | 268 ++++++++++++++++++++++ cartesian_grid_3d.py | 6 +- cmip5_all_plots.py | 45 ---- cmip5_atmos_climatology_netcdf.py | 182 --------------- cmip5_ecco2_rms_errors.py | 172 -------------- cmip5_eraint_rms_errors.py | 178 -------------- cmip5_field.py | 181 --------------- cmip5_max_uwind.py | 259 --------------------- cmip5_ocean_climatology_netcdf.py | 370 ------------------------------ cmip5_paths.py | 94 -------- ecco2_climatology_netcdf.py | 67 ------ ecco2_field.py | 71 ------ eraint_climatology_netcdf.py | 56 ----- eraint_field.py | 99 -------- file_guide.txt | 245 +++----------------- freezingpt_slice.py | 182 +++++++++++++++ i_slice.py | 181 +++++++++++++++ ismr_map.py | 8 +- massloss_map.py | 10 +- mmm_atmos_netcdf.py | 100 -------- mmm_ocean_netcdf.py | 94 -------- rignot_data | 25 -- 23 files changed, 685 insertions(+), 2214 deletions(-) create mode 100644 .gitignore create mode 100644 aice_hi_seasonal.py delete mode 100644 cmip5_all_plots.py delete mode 100644 cmip5_atmos_climatology_netcdf.py delete mode 100644 cmip5_ecco2_rms_errors.py delete mode 100644 cmip5_eraint_rms_errors.py delete mode 100644 cmip5_field.py delete mode 100644 cmip5_max_uwind.py delete mode 100644 cmip5_ocean_climatology_netcdf.py delete mode 100644 cmip5_paths.py delete mode 100644 ecco2_climatology_netcdf.py delete mode 100644 ecco2_field.py delete mode 100644 eraint_climatology_netcdf.py delete mode 100644 eraint_field.py create mode 100644 freezingpt_slice.py create mode 100644 i_slice.py delete mode 100644 mmm_atmos_netcdf.py delete mode 100644 mmm_ocean_netcdf.py delete mode 100644 rignot_data diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..d8976a7 --- /dev/null +++ b/.gitignore @@ -0,0 +1,6 @@ +*.pyc +*.png +*.log +*.jnl +area_comparison +rignot_data \ No newline at end of file diff --git a/aice_hi_seasonal.py b/aice_hi_seasonal.py new file mode 100644 index 0000000..96e7409 --- /dev/null +++ b/aice_hi_seasonal.py @@ -0,0 +1,268 @@ +from numpy import * +from netCDF4 import Dataset, num2date +from matplotlib.pyplot import * + +# Creates a 4x2 plot of seasonally averaged sea ice concentration (top row) and +# thickness (bottom row) over the last year of simulation. +# Input: +# cice_file = path to CICE output file with 5-day averages, containing at least +# one complete Dec-Nov period (if there are multiple such periods, +# this script uses the last one) +# 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 the figure +def aice_hi_seasonal (cice_file, save=False, fig_name=None): + + # Starting and ending months (1-based) for each season + start_month = [12, 3, 6, 9] + end_month = [2, 5, 8, 11] + # Starting and ending days of the month (1-based) for each season + # Assume no leap years, we'll fix this later if we need + start_day = [1, 1, 1, 1] + end_day = [28, 31, 31, 30] + # Number of days in each season (again, ignore leap years for now) + ndays_season = [90, 92, 92, 91] + # Season names for titles + season_names = ['DJF', 'MAM', 'JJA', 'SON'] + # Degrees to radians conversion + deg2rad = pi/180.0 + + # Read CICE grid and time values + id = Dataset(cice_file, 'r') + cice_lon = id.variables['TLON'][:-15,:] + cice_lat = id.variables['TLAT'][:-15,:] + time_id = id.variables['time'] + # Get the year, month, and day (all 1-based) for each output step + # These are 5-day averages marked with the last day's date. + cice_time = num2date(time_id[:], units=time_id.units, calendar=time_id.calendar.lower()) + + # Loop backwards through time indices to find the last one we care about + # (which contains 30 November in its averaging period) + end_t = -1 # Missing value flag + for t in range(size(cice_time)-1, -1, -1): + if cice_time[t].month == end_month[-1] and cice_time[t].day == end_day[-1]: + end_t = t + break + if cice_time[t].month == start_month[0] and cice_time[t].day in range(start_day[0], start_day[0]+4): + end_t = t + break + # Make sure we actually found it + if end_t == -1: + print 'Error: ' + cice_file + ' does not contain a complete Dec-Nov period' + return + + # Continue looping backwards to find the first time index we care about + # (which contains 1 December the previous year in its averaging period) + start_t = -1 # Missing value flag + for t in range(end_t-60, -1, -1): + if cice_time[t].month == start_month[0] and cice_time[t].day in range(start_day[0], start_day[0]+5): + start_t = t + break + # Make sure we actually found it + if start_t == -1: + print 'Error: ' + cice_file + ' does not contain a complete Dec-Nov period' + return + + # Check if end_t occurs on a leap year + leap_year = False + if mod(cice_time[end_t].year, 4) == 0: + # Years divisible by 4 are leap years + leap_year = True + if mod(cice_time[end_t].year, 100) == 0: + # Unless they're also divisible by 100, in which case they aren't + # leap years + leap_year = False + if mod(cice_time[end_t].year, 400) == 0: + # Unless they're also divisible by 400, in which case they are + # leap years after all + leap_year = True + if leap_year: + # Update last day in February + end_day[0] += 1 + ndays_season[0] += 1 + + # Initialise seasonal averages of CICE output + aice = ma.empty([4, size(cice_lon,0), size(cice_lon,1)]) + aice[:,:,:] = 0.0 + hi = ma.empty([4, size(cice_lon,0), size(cice_lon,1)]) + hi[:,:,:] = 0.0 + # Process one season at a time + for season in range(4): + season_days = 0 # Number of days in season; this will be incremented + next_season = mod(season+1, 4) + + # Find starting timestep + start_t_season = -1 + for t in range(start_t, end_t+1): + if cice_time[t].month == start_month[season] and cice_time[t].day in range(start_day[season], start_day[season]+5): + start_t_season = t + break + # Make sure we actually found it + if start_t_season == -1: + print 'Error: could not find starting timestep for season ' + season_names[season] + return + + # Find ending timestep + end_t_season = -1 + for t in range(start_t_season+1, end_t+1): + if cice_time[t].month == end_month[season] and cice_time[t].day == end_day[season]: + end_t_season = t + break + if cice_time[t].month == start_month[next_season] and cice_time[t].day in range(start_day[next_season], start_day[next_season]+4): + end_t_season = t + break + # Make sure we actually found it + if end_t_season == -1: + print 'Error: could not find ending timestep for season ' + season_names[season] + return + + # Figure out how many of the 5 days averaged in start_t_season are + # actually within this season + if cice_time[start_t_season].month == start_month[season] and cice_time[start_t_season].day == start_day[season] + 4: + # Starting day is in position 1 of 5; we care about all of them + start_days = 5 + elif cice_time[start_t_season].month == start_month[season] and cice_time[start_t_season].day == start_day[season] + 3: + # Starting day is in position 2 of 5; we care about the last 4 + start_days = 4 + elif cice_time[start_t_season].month == start_month[season] and cice_time[start_t_season].day == start_day[season]+ 2: + # Starting day is in position 3 of 5; we care about the last 3 + start_days = 3 + elif cice_time[start_t_season].month == start_month[season] and cice_time[start_t_season].day == start_day[season] + 1: + # Starting day is in position 4 of 5; we care about the last 2 + start_days = 2 + elif cice_time[start_t_season].month == start_month[season] and cice_time[start_t_season].day == start_day[season]: + # Starting day is in position 5 of 5; we care about the last 1 + start_days = 1 + else: + print 'Error for season ' + season_names[season] + ': starting index is month ' + str(cice_time[start_t_season].month) + ', day ' + str(cice_time[start_t_season].day) + return + + # Start accumulating data weighted by days + aice[season,:,:] += id.variables['aice'][start_t_season,:-15,:]*start_days + hi[season,:,:] += id.variables['hi'][start_t_season,:-15,:]*start_days + season_days += start_days + + # Beween start_t_season and end_t_season, we want all the days + for t in range(start_t_season+1, end_t_season): + aice[season,:,:] += id.variables['aice'][t,:-15,:]*5 + hi[season,:,:] += id.variables['hi'][t,:-15,:]*5 + season_days += 5 + + # Figure out how many of the 5 days averaged in end_t_season are + # actually within this season + if cice_time[end_t_season].month == start_month[next_season] and cice_time[end_t_season].day == start_day[next_season] + 3: + # Ending day is in position 1 of 5; we care about the first 1 + end_days = 1 + elif cice_time[end_t_season].month == start_month[next_season] and cice_time[end_t_season].day == start_day[next_season] + 2: + # Ending day is in position 2 of 5; we care about the first 2 + end_days = 2 + elif cice_time[end_t_season].month == start_month[next_season] and cice_time[end_t_season].day == start_day[next_season] + 1: + # Ending day is in position 3 of 5; we care about the first 3 + end_days = 3 + elif cice_time[end_t_season].month == start_month[next_season] and cice_time[end_t_season].day == start_day[next_season]: + # Ending day is in position 4 of 5; we care about the first 4 + end_days = 4 + elif cice_time[end_t_season].month == end_month[season] and cice_time[end_t_season].day == end_day[season]: + # Ending day is in position 5 of 5; we care about all 5 + end_days = 5 + else: + print 'Error for season ' + season_names[season] + ': ending index is month ' + str(cice_time[end_t_season].month) + ', day ' + str(cice_time[end_t_season].day) + return + + aice[season,:,:] += id.variables['aice'][end_t_season,:-15,:]*end_days + hi[season,:,:] += id.variables['hi'][end_t_season,:-15,:]*end_days + season_days += end_days + + # Check that we got the correct number of days + if season_days != ndays_season[season]: + print 'Error: found ' + str(season_days) + ' days instead of ' + str(ndays_season[season]) + return + + # Finished accumulating data, now convert from sum to average + aice[season,:,:] /= season_days + hi[season,:,:] /= season_days + + # Finished reading all CICE data + id.close() + + # Convert to spherical coordinates + cice_x = -(cice_lat+90)*cos(cice_lon*deg2rad+pi/2) + cice_y = (cice_lat+90)*sin(cice_lon*deg2rad+pi/2) + + # Set consistent colour levels + lev1 = linspace(0, 1, num=50) + lev2 = linspace(0, 3, num=50) + + # Plot + fig = figure(figsize=(20,9)) + # Loop over seasons + for season in range(4): + # aice + ax = fig.add_subplot(2, 4, season+1, aspect='equal') + img = contourf(cice_x, cice_y, aice[season,:,:], lev1, extend='both') + if season == 0: + text(-39, 0, 'aice (%)', fontsize=21, ha='right') + title(season_names[season], fontsize=24) + xlim([-35, 35]) + ylim([-33, 37]) + axis('off') + 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) + # hi + ax = fig.add_subplot(2, 4, season+5, aspect='equal') + img = contourf(cice_x, cice_y, hi[season,:,:], lev2, extend='both') + if season == 0: + text(-39, 0, 'hi (m)', fontsize=21, ha='right') + xlim([-35, 35]) + ylim([-33, 37]) + axis('off') + if season == 3: + cbaxes2 = fig.add_axes([0.92, 0.15, 0.01, 0.3]) + cbar2 = colorbar(img, ticks=arange(0,3+1,1), cax=cbaxes2) + cbar2.ax.tick_params(labelsize=16) + # Decrease space between plots + 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_hi_seasonal(cice_file, save, fig_name) + + + + + + + + + + + + + + + + + + + + + + diff --git a/cartesian_grid_3d.py b/cartesian_grid_3d.py index dc51e0c..ccc5bbe 100644 --- a/cartesian_grid_3d.py +++ b/cartesian_grid_3d.py @@ -31,8 +31,12 @@ def cartesian_grid_3d (lon, lat, h, zice, theta_s, theta_b, hc, N, zeta=None): # bottom edges of each cell z_edges = zeros((N+1, num_lat, num_lon)) z_edges[1:-1,:,:] = 0.5*(z[0:-1,:,:] + z[1:,:,:]) - # At surface, z=zice; at bottom, extrapolate + # At surface, z=zice z_edges[-1,:,:] = zice[:,:] + # Add zeta if it exists + if zeta is not None: + z_edges[-1,:,:] += zeta[:,:] + # At bottom, extrapolate z_edges[0,:,:] = 2*z[0,:,:] - z_edges[1,:,:] # Now find dz dz = z_edges[1:,:,:] - z_edges[0:-1,:,:] diff --git a/cmip5_all_plots.py b/cmip5_all_plots.py deleted file mode 100644 index 657b024..0000000 --- a/cmip5_all_plots.py +++ /dev/null @@ -1,45 +0,0 @@ -from cmip5_paths import * -from cmip5_plot import * - -# Call cmip5_plot.py for all models (including the multi-model mean), all -# variables (atmosphere and ocean), and all seasons. -# Note that in order to run this script, you must first run: -# eraint_climatology_netcdf.py, ecco2_climatology_netcdf.py, -# cmip5_atmos_climatology_netcdf.py, cmip5_ocean_climatology_netcdf.py, -# mmm_atmos_netcdf.py, and mmm_ocean_netcdf.py -# to generate the necessary NetCDF files. Also, the directory "cmip5/" must -# exist. -# Input: season = string containing the season key: 'djf', 'mam', 'jja', 'son', -# or 'annual' -def cmip5_all_plots (season): - - # All possible variable names - var_names = ['Pair', 'Tair', 'Hair', 'cloud', 'Uwind', 'Vwind', 'precip', 'snow', 'evap', 'swrad', 'lwrad', 'temp', 'salt', 'u', 'v'] - - # Make a list of all possible model names - models = build_model_list() - model_names = [] - for model in models: - model_names.append(model.name) - model_names.append('MMM') - - # Call cmip5_plot for each variable - for i in range(len(var_names)): - var = var_names[i] - print 'Plotting ' + var - cmip5_plot(var, season, model_names, True, 'cmip5/' + var + '_' + season + '.png') - - -# Command-line interface -if __name__ == "__main__": - - # Run for each season key at a time - cmip5_all_plots('djf') - cmip5_all_plots('mam') - cmip5_all_plots('jja') - cmip5_all_plots('son') - cmip5_all_plots('annual') - - - - diff --git a/cmip5_atmos_climatology_netcdf.py b/cmip5_atmos_climatology_netcdf.py deleted file mode 100644 index 5f3ea6e..0000000 --- a/cmip5_atmos_climatology_netcdf.py +++ /dev/null @@ -1,182 +0,0 @@ -from cmip5_paths import * -from cmip5_field import * -from eraint_field import * -from numpy import * -from netCDF4 import Dataset -from scipy.interpolate import RegularGridInterpolator - -# NB for raijin users: RegularGridInterpolator needs python/2.7.6 but the -# default is 2.7.3. Before running this script, switch them as follows: -# module unload python/2.7.3 -# module unload python/2.7.3-matplotlib -# module load python/2.7.6 -# module load python/2.7.6-matplotlib - -# For the given CMIP5 model, calculates the monthly climatology from 1992-2005 -# inclusive for 11 atmospheric variables. Interpolates to the ERA-Interim grid -# and saves to a NetCDF file. Note that to run this script, you must have -# previously run eraint_climatology_netcdf.py. -# Input: -# model_name = name of model (this must match the list in cmip5_paths.py) -def cmip5_atmos_climatology_netcdf (model_name): - - # Experiment name - expt = 'historical' - # Years over which to calculate climatology - start_year = 1992 - end_year = 2005 - # CMIP5 variable names - var_names = ['ps', 'tas', 'huss', 'clt', 'uas', 'vas', 'pr', 'prsn', 'evspsbl', 'rsds', 'rlds'] - # Variable names to use in NetCDF file - var_names_output = ['Pair', 'Tair', 'Hair', 'cloud', 'Uwind', 'Vwind', 'precip', 'snow', 'evap', 'swrad', 'lwrad'] - # Units of final variables (note some conversions in cmip5_field) - var_units = ['kPa', 'degC', '1', '%', 'm/s', 'm/s', '10^6 kg/m^2/s', '10^6 kg/m^2/s', '10^6 kg/m^2/s', 'W/m^2', 'W/m^2'] - # Path to output NetCDF file - output_file = '/short/y99/kaa561/CMIP5_forcing/atmos/' + model_name + '.nc' - # Path to corresponding ERA-Interim file (created using - # eraint_climatology_netcdf.py) - eraint_file = '/short/y99/kaa561/CMIP5_forcing/atmos/ERA-Interim.nc' - - # Read ERA-Interim grid - id = Dataset(eraint_file, 'r') - era_lon = id.variables['longitude'][:] - era_lat = id.variables['latitude'][:] - id.close() - - # Build a list of Model objects for CMIP5 - all_models = build_model_list() - # Build a corresponding list of all CMIP5 model names - all_model_names = [] - for model in all_models: - all_model_names.append(model.name) - # Find index of model_name in all_model_names, and select the Model object - # at the same index of all_models - # Now model_name has a corresponding Model object - model = all_models[all_model_names.index(model_name)] - - # Loop over variables - for i in range(len(var_names)): - var = var_names[i] - print 'Processing variable ' + var - - # Read monthly climatology for this variable - model_data, model_lon, model_lat, tmp = cmip5_field(model, expt, var, start_year, end_year) - # Set up array for climatology interpolated to ERA-Interim grid - model_data_interp = ma.empty([12, size(era_lat), size(era_lon)]) - if model_data is not None: - # Interpolate one month at a time - for t in range(size(model_data,0)): - model_data_interp[t,:,:] = interp_model2era(model_data[t,:,:], model_lon, model_lat, era_lon, era_lat) - else: - # No data (missing variable in CMIP5 archive) - model_data_interp[:,:,:] = ma.masked - if i == 0: - # Set up NetCDF file on the first iteration - print 'Setting up ' + output_file - id = Dataset(output_file, 'w') - # Define dimensions - id.createDimension('longitude', size(era_lon)) - id.createDimension('latitude', size(era_lat)) - id.createDimension('time', 12) - # Define dimension variables and fill with data - id.createVariable('longitude', 'f8', ('longitude')) - id.variables['longitude'].units = 'degrees' - id.variables['longitude'][:] = era_lon - id.createVariable('latitude', 'f8', ('latitude')) - id.variables['latitude'].units = 'degrees' - id.variables['latitude'][:] = era_lat - id.createVariable('time', 'f8', ('time')) - id.variables['time'].units = 'month' - id.variables['time'][:] = arange(1, 12+1) - # Define new variable and fill with data - id.createVariable(var_names_output[i], 'f8', ('time', 'latitude', 'longitude')) - id.variables[var_names_output[i]].units = var_units[i] - id.variables[var_names_output[i]][:,:,:] = model_data_interp - id.close() - - -# Interpolate a given CMIP5 field to the ERA-Interim grid. -# Input: -# model_data = 2D array (size mxn) of model data for the given variable -# model_lon = 1D array (length n) of longitude for this model -# model_lat = 1D array (length m) of latitude for this model -# era_lon = 1D array (length q) of longitude on the ERA-Interim grid -# era_lat = 1D array (length p) of latitude on the ERA-Interim grid -# Output: -# data_interp = 2D array (size pxq) of model data interpolated to the -# ERA-Interim grid -def interp_model2era(model_data, model_lon, model_lat, era_lon, era_lat): - - # Make sure the model's longitude goes from 0 to 360, not -180 to 180 - index = model_lon < 0 - model_lon[index] = model_lon[index] + 360 - - # CMIP5 model axes don't wrap around; there is a gap between almost-180W - # and almost-180E (these are 0 to 360 but you get the idea) and depending - # on the grid, we may need to interpolate in this gap. - # So copy the last longitude value (mod 360) to the beginning, and the - # first longitude value (mod 360) to the end. - model_lon_wrap = zeros(size(model_lon)+2) - model_lon_wrap[0] = model_lon[-1] - 360 - model_lon_wrap[1:-1] = model_lon - model_lon_wrap[-1] = model_lon[0] + 360 - model_lon = model_lon_wrap - # Copy the westernmost and easternmost data points to match - model_data_wrap = ma.array(zeros((size(model_lat), size(model_lon)))) - model_data_wrap[:,1:-1] = model_data - model_data_wrap[:,0] = model_data[:,-1] - model_data_wrap[:,-1] = model_data[:,0] - model_data = model_data_wrap - - if amin(model_lat) > amin(era_lat): - # Add a point at 90S - model_lat_new = zeros(size(model_lat)+1) - model_data_new = ma.array(zeros((size(model_lat_new), size(model_lon)))) - if model_lat[0] > model_lat[1]: - model_lat_new[0:-1] = model_lat - model_lat_new[-1] = -90.0 - model_data_new[0:-1,:] = model_data - model_data_new[-1,:] = model_data[-1,:] - elif model_lat[0] < model_lat[1]: - model_lat_new[1:] = model_lat - model_lat_new[0] = -90.0 - model_data_new[1:,:] = model_data - model_data_new[0,:] = model_data[0,:] - model_lat = model_lat_new - model_data = model_data_new - if amax(model_lat) < amax(era_lat): - # Add a point at 90N - model_lat_new = zeros(size(model_lat)+1) - model_data_new = ma.array(zeros((size(model_lat_new), size(model_lon)))) - if model_lat[0] > model_lat[1]: - model_lat_new[1:] = model_lat - model_lat_new[0] = 90.0 - model_data_new[1:,:] = model_data - model_data_new[0,:] = model_data[0,:] - elif model_lat[0] < model_lat[1]: - model_lat_new[0:-1] = model_lat - model_lat_new[-1] = 90.0 - model_data_new[0:-1,:] = model_data - model_data_new[-1,:] = model_data[-1,:] - model_lat = model_lat_new - model_data = model_data_new - - # Get 2D mesh of ERA-Interim lat and lon - era_lon_2d, era_lat_2d = meshgrid(era_lon, era_lat) - # Build an interpolation function for model_data - interp_function = RegularGridInterpolator((model_lat, model_lon), model_data) - # Call it for the ERA-Interim grid - data_interp = interp_function((era_lat_2d, era_lon_2d)) - - return data_interp - - -# Command-line interface -if __name__ == "__main__": - - # Process one model at a time - models = build_model_list() - for model in models: - print model.name - cmip5_atmos_climatology_netcdf(model.name) - diff --git a/cmip5_ecco2_rms_errors.py b/cmip5_ecco2_rms_errors.py deleted file mode 100644 index 7974fbe..0000000 --- a/cmip5_ecco2_rms_errors.py +++ /dev/null @@ -1,172 +0,0 @@ -from cmip5_paths import * -from numpy import * -from netCDF4 import Dataset -from matplotlib.pyplot import * - -# Calculate root-mean square errors (adapted from Gleckler et al., 2008 to be -# in depth-longitude space instead of latitude-longitude space) for each of 39 -# CMIP5 models and the multi-model mean, with respect to 4 ECCO2 ocean -# variables at the northern boundary of ROMS (currently 30S). Use the monthly -# climatology averaged over 1992-2005 inclusive. Also calculate the relative -# errors as in Gleckler et al. and make a "portrait plot" of coloured tiles -# in a model vs variable matrix. Save both rms errors and relative errors into -# text files. Note that to run this script, you must have previously run -# ecco2_climatology_netcdf.py, cmip5_ocean_climatology_netcdf.py, and -# mmm_ocean_netcdf.py. -def cmip5_ecco2_rms_errors (): - - # Path to directory containing climatolgoy NetCDF files (created by - # ecco2_climatology_netcdf.py and cmip5_ocean_climatology_netcdf.py) - directory = '/short/y99/kaa561/CMIP5_forcing/ocean/' - # Variable names in NetCDF files - var_names = ['temp', 'salt', 'u', 'v'] - # Radius of the Earth in metres - r = 6.371e6 - # Degrees to radians conversion factor - deg2rad = pi/180.0 - # Latitude of the northern boundary - nbdry = -30.0 - - # Read ECCO2 grid - id = Dataset(directory + 'ECCO2.nc', 'r') - lon = id.variables['longitude'][:] - depth = id.variables['depth'][:] - id.close() - - # Calculate dx, dz, dt - # ECCO2 has constant spacing for longitude; find this value (should be 0.25 - # degrees) - dlon = mean(lon[1:] - lon[0:-1]) - # dx = r*cos(lat)*dlon where lat and dlon are converted to radians - dx = abs(r*cos(nbdry*deg2rad)*dlon*deg2rad) - # We have depth values in the middle of each cell; interpolate to the edges - z_edges = zeros(size(depth)+1) - z_edges[1:-1] = 0.5*(depth[0:-1] + depth[1:]) - # Assume surface is 0 - z_edges[0] = 0.0 - # Extrapolate to bottom - z_edges[-1] = 2*depth[-1] - z_edges[-2] - # Now calculate dz - dz = z_edges[1:] - z_edges[0:-1] - # dt = number of seconds in each month - dt = array([31.0, 28.25, 31.0, 30.0, 31.0, 30.0, 31.0, 31.0, 30.0, 31.0, 30.0, 31.0])*24*60*60 - # Copy dz and dt to 3D arrays with dimension time x depth x longitude - # We don't have to do this for dx because it's a scalar - dz = transpose(tile(dz, (size(lon), 1))) - dz = tile(dz, (12, 1, 1)) - dt = tile(dt, (size(lon), size(depth), 1)) - dt = transpose(dt) - - # Build a list of CMIP5 Model objects - models = build_model_list() - # Build a corresponding list of model names - model_names = [] - for model in models: - model_names.append(model.name) - # Add the multi-model mean to the list - model_names.append('MMM') - - # Set up arrays for RMS errors and relative errors, of dimension - # variable x model - rms_errors = ma.empty([len(var_names), len(model_names)]) - relative_errors = ma.empty(shape(rms_errors)) - - # Loop over variables - for i in range(len(var_names)): - var = var_names[i] - print 'Variable ' + var - - # Read ECCO2 data - print 'Processing ECCO2' - id = Dataset(directory + 'ECCO2.nc', 'r') - ecco_data = id.variables[var][:,:,:] - id.close() - - if i == 0: - # Mask dz with the land mask on the first iteration - # If at least one of the integrands are masked, land points - # will be automatically excluded from the integrals - dz_mask = ma.masked_where(ecco_data.mask, dz) - - # Loop over models - for j in range(len(model_names)): - model_name = model_names[j] - print 'Processing ' + model_name - - # Read model data - # This has already been interpolated to the ECCO2 grid at the - # northern boundary of ROMS - id = Dataset(directory + model_name + '.nc', 'r') - model_data = id.variables[var][:,:,:] - id.close() - - if all(model_data.mask): - # This variable is missing; mask out the RMS error - rms_errors[i,j] = ma.masked - else: - # Calculate RMS error: square root of the sum of the squares of - # the difference between this model and ECCO2 at every point, - # weighted by dx*dz*dt, with land points excluded since dz - # has been masked - rms_errors[i,j] = sqrt(sum((model_data - ecco_data)**2*dx*dz_mask*dt)/sum(dx*dz_mask*dt)) - - # Calculate relative error as (E-E')/E where E is the RMS error for each - # model, and E' is the median error across all models for this variable - median_error = median(rms_errors[i,:-1]) - for j in range(len(models)+1): - relative_errors[i,j] = (rms_errors[i,j] - median_error)/rms_errors[i,j] - - # Make the portrait plot - fig = figure(figsize=(8,20)) - ax = fig.add_subplot(111) - img = ax.pcolormesh(transpose(relative_errors), vmin=-1.5, vmax=1.5, cmap='RdBu_r') - - # Configure plot - # Add model names and variable names to the axes - xticks(arange(0.5, len(var_names)+0.5), var_names) - yticks(arange(0.5, len(model_names)+0.5), model_names) - xlim([0, len(var_names)]) - ylim([len(model_names), 0]) - title('Relative error', fontsize=16) - # Move the plot over a bit to make room for all the model names - box = ax.get_position() - ax.set_position([box.x0+box.width*0.1, box.y0, box.width*0.9, box.height]) - cbaxes = fig.add_axes([0.3,0.03,0.5,0.025]) - colorbar(img, orientation='horizontal', ticks=arange(-1.5, 1.75, 0.5), cax=cbaxes) - savefig('relative_errors_ecco2.png') - - # Write RMS errors to a text file, tabulated nicely - f = open('rms_errors_ecco2.txt', 'w') - f.write('{0:16}'.format('Model')) - for var in var_names: - f.write('{0:16}'.format(var)) - f.write('\n') - for j in range(len(model_names)): - f.write('{0:16}'.format(model_names[j])) - for i in range(len(var_names)): - f.write('{0:16}'.format(str(round(rms_errors[i,j],4)))) - f.write('\n') - f.close() - - # Write relative errors to a text file, tabulated nicely - f = open('relative_errors_ecco2.txt', 'w') - f.write('{0:16}'.format('Model')) - for var in var_names: - f.write('{0:16}'.format(var)) - f.write('\n') - for j in range(len(model_names)): - f.write('{0:16}'.format(model_names[j])) - for i in range(len(var_names)): - f.write('{0:16}'.format(str(round(relative_errors[i,j],4)))) - f.write('\n') - f.close() - - -# Command-line interface -if __name__ == '__main__': - - cmip5_ecco2_rms_errors() - - - - diff --git a/cmip5_eraint_rms_errors.py b/cmip5_eraint_rms_errors.py deleted file mode 100644 index bd60617..0000000 --- a/cmip5_eraint_rms_errors.py +++ /dev/null @@ -1,178 +0,0 @@ -from cmip5_paths import * -from numpy import * -from netCDF4 import Dataset -from matplotlib.pyplot import * - -# Calculate root-mean-square errors (as in Gleckler et al., 2008) for each of -# 39 CMIP5 models and the multi-model mean, with respect to 11 ERA-Interim -# atmospheric variables. The domain is the Southern Ocean (all longitudes, and -# latitudes from the northern boundary of ROMS to the southernmost ocean point -# not in an ice shelf cavity) and the monthly climatology averaged over -# 1992-2005 inclusive. Also calculate the relative errors as in Gleckler et al. -# and make a "portrait plot" of coloured tiles in a model vs variable matrix. -# Save both rms errors and relative errors into text files. -# Note that to run this script, you must have previously run -# eraint_climatology_netcdf.py, cmip5_atmos_climatology_netcdf.py, and -# mmm_atmos_netcdf.py. -def cmip5_eraint_rms_errors(): - - # Path to directory containing climatology NetCDF files (created by - # eraint_climatology_netcdf.py and cmip5_atmos_climatology_netcdf.py) - directory = '/short/y99/kaa561/CMIP5_forcing/atmos/' - # Variable names in NetCDF files - var_names = ['Pair', 'Tair', 'Hair', 'cloud', 'Uwind', 'Vwind', 'precip', 'snow', 'evap', 'swrad', 'lwrad'] - # Path to ROMS grid - roms_grid = '/short/m68/kaa561/ROMS-CICE-MCT/apps/common/grid/circ30S_quarterdegree_10m.nc' - # Radius of the Earth in metres - r = 6.371e6 - # Degrees to radians conversion factor - deg2rad = pi/180.0 - - # Read ROMS grid - id = Dataset(roms_grid, 'r') - lat_rho = id.variables['lat_rho'][:,:] - mask_rho = id.variables['mask_rho'][:,:] - mask_zice = id.variables['mask_zice'][:,:] - id.close() - # Determine latitude bounds of ocean points not in ice shelf cavities - lat_roms = ma.masked_where(mask_rho-mask_zice==0, lat_rho) - latS = amin(lat_roms) - latN = amax(lat_roms) - - # Read ERA-Interim grid - id = Dataset(directory + 'ERA-Interim.nc', 'r') - lon = id.variables['longitude'][:] - lat = id.variables['latitude'][:] - id.close() - - # Find bounds on j to extend just past latS and latN - j_min = nonzero(lat < latN)[0][0] - 1 - j_max = nonzero(lat < latS)[0][0] + 1 - lat = lat[j_min:j_max] - - # Calculate dx, dy, dt - # ERA-Interim has constant spacing for both lat and lon; find these values - # (should be 0.75 degrees) - dlon = mean(lon[1:] - lon[0:-1]) - dlat = mean(lat[1:] - lat[0:-1]) - # dx = r*cos(lat)*dlon where lat and dlon are converted to radians - dx = abs(r*cos(lat*deg2rad)*dlon*deg2rad) - # dy = r*dlat where dlat is converted to radians - dy = abs(r*dlat*deg2rad) - # dt = number of seconds in each month - dt = array([31.0, 28.25, 31.0, 30.0, 31.0, 30.0, 31.0, 31.0, 30.0, 31.0, 30.0, 31.0])*24*60*60 - # Copy dx and dt to 3D arrays with dimension time x latitude x longitude - # We don't have to do this for dy because it's a scalar - dx = transpose(tile(dx, (size(lon), 1))) - dx = tile(dx, (12, 1, 1)) - dt = tile(dt, (size(lon), size(lat), 1)) - dt = transpose(dt) - - # Build a list of CMIP5 Model objects - models = build_model_list() - # Build a corresponding list of model names - model_names = [] - for model in models: - model_names.append(model.name) - # Add the multi-model mean to the list - model_names.append('MMM') - - # Set up arrays for RMS errors and relative errors, of dimension - # variable x model - rms_errors = ma.empty([len(var_names), len(model_names)]) - relative_errors = ma.empty(shape(rms_errors)) - - # Loop over variables - for i in range(len(var_names)): - var = var_names[i] - print 'Variable ' + var - - # Read ERA-Interim data, trimmed to latitude bounds - print 'Processing ERA-Interim' - id = Dataset(directory + 'ERA-Interim.nc', 'r') - era_data = id.variables[var][:,j_min:j_max,:] - id.close() - - # Loop over models - for j in range(len(model_names)): - model_name = model_names[j] - print 'Processing ' + model_name - - # Read model data, trimmed to latitude bounds - # This has already been interpolated to the ERA-Interim grid - id = Dataset(directory + model_name + '.nc', 'r') - model_data = id.variables[var][:,j_min:j_max,:] - id.close() - - # Figure out if this variable is missing - try: - mask = model_data.mask - except(AttributeError): - # There is no mask at all - mask = False - if all(mask): - # Everything is masked, i.e. the variable is missing - # Mask out the RMS error too - rms_errors[i,j] = ma.masked - else: - # Calculate RMS error: square root of the sum of the squares - # of the difference between this model and ERA-Interim at every - # point, weighted by dx*dy*dt - rms_errors[i,j] = sqrt(sum((model_data - era_data)**2*dx*dy*dt)/sum(dx*dy*dt)) - - # Calculate relative error as (E-E')/E where E is the RMS error for each - # model, and E' is the median error across all models for this variable - median_error = median(rms_errors[i,:-1]) - for j in range(len(models)+1): - relative_errors[i,j] = (rms_errors[i,j] - median_error)/rms_errors[i,j] - - # Make the portrait plot - fig = figure(figsize=(8,20)) - ax = fig.add_subplot(111) - img = ax.pcolormesh(transpose(relative_errors), vmin=-0.8, vmax=0.8, cmap='RdBu_r') - - # Configure plot - # Add model names and variable names to the axes - xticks(arange(0.5, len(var_names)+0.5), var_names, rotation=-45) - yticks(arange(0.5, len(model_names)+0.5), model_names) - xlim([0, len(var_names)]) - ylim([len(model_names), 0]) - title('Relative error', fontsize=16) - # Move the plot over a bit to make room for all the model names - box = ax.get_position() - ax.set_position([box.x0+box.width*0.1, box.y0, box.width*0.9, box.height]) - cbaxes = fig.add_axes([0.3,0.03,0.5,0.025]) - colorbar(img, orientation='horizontal', ticks=arange(-0.75, 1, 0.25), cax=cbaxes) - savefig('relative_errors_eraint.png') - - # Write RMS errors to a text file, tabulated nicely - f = open('rms_errors_eraint.txt', 'w') - f.write('{0:16}'.format('Model')) - for var in var_names: - f.write('{0:16}'.format(var)) - f.write('\n') - for j in range(len(model_names)): - f.write('{0:16}'.format(model_names[j])) - for i in range(len(var_names)): - f.write('{0:16}'.format(str(round(rms_errors[i,j],4)))) - f.write('\n') - f.close() - - # Write relative errors to a text file, tabulated nicely - f = open('relative_errors_eraint.txt', 'w') - f.write('{0:16}'.format('Model')) - for var in var_names: - f.write('{0:16}'.format(var)) - f.write('\n') - for j in range(len(model_names)): - f.write('{0:16}'.format(model_names[j])) - for i in range(len(var_names)): - f.write('{0:16}'.format(str(round(relative_errors[i,j],4)))) - f.write('\n') - f.close() - - -# Command-line interface -if __name__ == '__main__': - - cmip5_eraint_rms_errors() diff --git a/cmip5_field.py b/cmip5_field.py deleted file mode 100644 index 179384b..0000000 --- a/cmip5_field.py +++ /dev/null @@ -1,181 +0,0 @@ -from numpy import * -from netCDF4 import Dataset, num2date, date2num -from os import listdir -from scipy.interpolate import interp1d - -# Read CMIP5 output for the given model, experiment, and variable name. Return -# the monthly climatology as well as the grid. -# Input: -# model = Model object (class definition in cmip5_paths.py) -# expt = string containing name of experiment, eg 'historical' -# var_name = string containing name of variable, eg 'tas' -# start_year, end_year = integers containing years over which to calculate -# monthly climatology, from the beginning of start_year -# to the end of end_year. Therefore, if -# start_year = end_year, this script will read one year -# of output with no climatological averaging. -# Output: -# data = array of model output, with dimension month x latitude x longitude -# (for atmosphere variables, at the surface) or month x depth x latitude -# x longitude (for ocean variables), possibly with units converted to be -# more comparable to other models and/or reanalyses -# lon, lat = longitude and latitude arrays (1D or 2D) -# depth = depth array (1D or 3D) if ocean variable, otherwise None -def cmip5_field (model, expt, var_name, start_year, end_year): - - # Conversion from K to C - degKtoC = -273.15 - - # Figure out whether it is an atmosphere or ocean variable - if var_name in ['ps', 'tas', 'huss', 'clt', 'uas', 'vas', 'pr', 'prsn', 'evspsbl', 'rsds', 'rlds']: - realm = 'atmos' - elif var_name in ['thetao', 'so', 'uo', 'vo', 'zos']: - realm = 'ocean' - else: - print 'Unknown variable' - # Exit early - return None, None, None, None - - # There is something weird going on with evaporation in the Norwegian GCMs; - # they claim to have units of kg/m^2/s but there's no way that's correct. - # Exclude them for now, until we figure out what the units actually are. - if var_name == 'evspsbl' and model.name in ['NorESM1-M', 'NorESM1-ME']: - print 'Skipping ' + model.name + ' because evaporation units are screwy' - # Exit early - return None, None, None, None - - # Get the directory where the model output is stored - path = model.get_directory(expt, var_name) - # If the string is empty, this output doesn't exist - if path == '': - print 'Warning: no data found for model ' + model.name + ', experiment ' + expt + ', variable ' + var_name - # Exit early - return None, None, None, None - - data = None - # Number of records for each month - num_months = zeros(12) - - # Loop over all files in this directory - for file in listdir(path): - - # Check every netCDF file - if file.endswith('.nc'): - - # Read the time values - id = Dataset(path + file, 'r') - time_id = id.variables['time'] - if amin(time_id[:]) < 0: - # Missing values here; this occurs for one 1900-1949 file - # We can just skip it - break - # Convert to datetime objects - curr_units = time_id.units - if curr_units == 'days since 0001-01': - curr_units = 'days since 0001-01-01' - if curr_units == 'days since 0000-01-01 00:00:00': - curr_units = 'days since 0001-01-01 00:00:00' - time = num2date(time_id[:], units=curr_units, calendar=time_id.calendar) - - # Check if the time values in this file actually contain any - # dates we're interested in - if time[0].year > end_year or time[-1].year < start_year: - # No overlap, skip this file and go to the next one - id.close() - else: - # Read the grid - lat = id.variables['lat'][:] - lon = id.variables['lon'][:] - if realm == 'ocean': - if model.name == 'inmcm4': - # inmcm4 has sigma coordinates; convert to z - sigma = id.variables['lev'][:] - h = id.variables['depth'][:,:] - depth = ma.empty([size(sigma), size(h,0), size(h,1)]) - for k in range(size(sigma)): - depth[k,:,:] = -sigma[k]*h - else: - depth = id.variables['lev'][:] - else: - depth = None - - # Read model output, one timestep at a time to save memory - for t in range(size(time)): - if realm == 'atmos': - data_tmp = id.variables[var_name][t,:,:] - elif realm == 'ocean': - data_tmp = id.variables[var_name][t,:,:,:] - # Some of the CMIP5 ocean models are not saved as masked - # arrays, but rather as regular arrays with the value 0 - # at all land points. Catch these with a try-except block - # and convert to masked arrays. - try: - mask = data_tmp.mask - except (AttributeError): - data_tmp = ma.masked_where(data_tmp==0, data_tmp) - - if data is None: - # Set up data array of correct size - if realm == 'atmos': - data = ma.empty([12, size(data_tmp,0), size(data_tmp,1)]) - for tt in range(12): - # Initialise as zeros masked with land mask - data[tt,:,:] = data_tmp[:,:]*0.0 - elif realm == 'ocean': - data = ma.empty([12, size(data_tmp,0), size(data_tmp,1), size(data_tmp,2)]) - for tt in range(12): - data[tt,:,:,:] = data_tmp[:,:,:]*0.0 - - # Figure out if it falls between the dates we want - if time[t].year >= start_year and time[t].year <= end_year: - # Find the month index, increment curr_month, and add - # to data array - curr_month = time[t].month-1 - num_months[curr_month] += 1 - if realm == 'atmos': - data[curr_month,:,:] += data_tmp[:,:] - elif realm == 'ocean': - data[curr_month,:,:,:] += data_tmp[:,:,:] - - # Check if we actually read any data - if data is None: - print 'No files found in specified date range' - # Exit early - return None, None, None, None - - # Convert from monthly sums to monthly averages - for t in range(12): - if num_months[t] == 0: - # None for this month; mask it - if realm == 'atmos': - data[t,:,:] = ma.masked - elif realm == 'ocean': - data[t,:,:,:] = ma.masked - else: - if realm == 'atmos': - data[t,:,:] /= num_months[t] - elif realm == 'ocean': - data[t,:,:,:] /= num_months[t] - - # Conversions if necessary - if var_name in ['pr', 'prsn', 'evspsbl']: - # Convert precip/snowfall/evap from kg/m^2/s to - # 1e-6 kg/m^2/s - data = 1e6*data - elif var_name == 'ps': - # Convert surface pressure from Pa to kPa - data = 1e-3*data - elif var_name == 'tas': - # Convert temperature from K to C - data = data + degKtoC - elif var_name == 'thetao' and amin(data) > 100: - # Convert ocean temperature from K to C if needed - data = data + degKtoC - elif var_name == 'so' and amax(data) < 1: - # Convert salinity from fraction to psu if needed - data = 1e3*data - - return data, lon, lat, depth - - - diff --git a/cmip5_max_uwind.py b/cmip5_max_uwind.py deleted file mode 100644 index ac26429..0000000 --- a/cmip5_max_uwind.py +++ /dev/null @@ -1,259 +0,0 @@ -from cmip5_paths import * -from numpy import * -from netCDF4 import Dataset -from matplotlib.pyplot import * -from matplotlib.font_manager import FontProperties - -# Create two plots: (1) the maximum zonal wind speed between 30S and 65S, and -# (2) the latitude of that maximum zonal wind speed, both against longitude. -# Plot results for ERA-Interim as well as the given CMIP5 models, averaged over -# the given season, for the 1992-2005 climatology. Note that in order to run -# this script, you must first run eraint_climatology_netcdf.py, -# cmip5_atmos_climatology_netcdf.py, and mmm_atmos_netcdf.py. -# Input: -# model_names = list of strings containing model names (either one of the CMIP5 -# models from cmip5_paths.py, or 'MMM' for the multi-model mean) -# season = string containing the season key: 'djf', 'mam', 'jja', 'son', or -# 'annual' -# save = optional boolean indicating that the plot should be saved to a file -# rather than displayed on the screen -# fig_names = if save=True, an array of length 2 containing the names of the -# files to save the 2 figures to -def cmip5_max_uwind (model_names, season, save=False, fig_names=None): - - # Directory where the NetCDF climatology files, creating using the scripts - # listed above, are stored - directory = '/short/y99/kaa561/CMIP5_forcing/atmos/' - # Bounds on latitude to search for the maximum uwind between - latS = -65 - latN = -30 - - # 40 visually distinct (or as visually distinct as possible) colours for - # plotting, generated using http://phrogz.net/css/distinct-colors.html - all_cmip5_colours = [(1,0,0), (0.58,0,0), (1,0.37,0.37), (0.58,0.22,0.22), (0.37,0.14,0.14), (1,0.74,0.74), (0.58,0.43,0.43), (1,0.78,0), (0.37,0.29,0), (0.58,0.5,0.22), (1,0.95,0.74), (0.58,0.55,0.43), (0.43,1,0), (0.16,0.37,0), (0.51,0.79,0.29), (0,0.58,0.2), (0.37,1,0.59), (0.74,1,0.83), (0.43,0.58,0.48), (0.27,0.37,0.31), (0.29,0.73,0.79), (0.14,0.34,0.37), (0.74,0.96,1), (0,0.08,1), (0,0.05,0.58), (0,0.03,0.37), (0.37,0.42,1), (0.22,0.24,0.58), (0.58,0.6,0.79), (0.55,0,0.79), (0.81,0.37,1), (0.47,0.22,0.58), (0.3,0.14,0.37), (0.92,0.74,1), (0.53,0.43,0.58), (0.34,0.27,0.37), (0.58,0,0.3), (1,0.37,0.69), (0.37,0.14,0.26), (1,0.74,0.87)] - # Figure out how many colours we need, and select them evenly from the list - num_colours = len(all_cmip5_colours) - num_models = len(model_names) - cmip5_colours = [] - for i in range(num_models): - index = int(ceil(i*float(num_colours)/num_models)) - cmip5_colours.append(all_cmip5_colours[index]) - - # Choose the month indices we are interested in, based on the season - if season == 'djf': - months = [12, 1, 2] - elif season == 'mam': - months = [3, 4, 5] - elif season == 'jja': - months = [6, 7, 8] - elif season == 'son': - months = [9, 10, 11] - elif season == 'annual': - months = range(1,12+1) - - # Initialise figures - fig1, ax1 = subplots(figsize=(12,9)) - fig2, ax2 = subplots(figsize=(12,9)) - - print 'Processing ERA-Interim' - # Read ERA-Interim grid - id = Dataset(directory + 'ERA-Interim.nc', 'r') - lon = id.variables['longitude'][:] - lat = id.variables['latitude'][:] - - # Find bounds on j to extend just past latS and latN - j_min = nonzero(lat < latN)[0][0] - 1 - j_max = nonzero(lat < latS)[0][0] + 1 - lat = lat[j_min:j_max] - - # Read ERA-Interim data - era_data = id.variables['Uwind'][:,j_min:j_max,:] - id.close() - for month in range(12): - # Mask out the months we don't care about - if month+1 not in months: - era_data[month,:,:] = ma.masked - # Take time average (this will automatically exclude the masked months) - era_data = mean(era_data, axis=0) - # Find the maximum along the latitude axis for each longitude point - era_max = amax(era_data, axis=0) - # Find the latitude indices of these maximum values - era_index = argmax(era_data, axis=0) - # Find the latitude values at these indices - era_loc = zeros(size(era_index)) - for i in range(size(era_index)): - era_loc[i] = lat[era_index[i]] - - # Add to plots (setting zorder means it will be on top) - ax1.plot(lon, era_max, label='ERA-Interim', color='black', linewidth=3, zorder=num_models+1) - ax2.plot(lon, era_loc, label='ERA-Interim', color='black', linewidth=3, zorder=num_models+1) - - # Loop over models - for j in range(len(model_names)): - model_name = model_names[j] - print 'Processing ' + model_name - # Read model data (note it has already been interpolated to ERA-Interim - # grid) - id = Dataset(directory + model_name + '.nc', 'r') - model_data = id.variables['Uwind'][:,j_min:j_max,:] - id.close() - # Check for missing data - try: - mask = model_data.mask - except(AttributeError): - # There is no mask - mask = False - if all(mask): - # Everything is masked, so the variable is missing - pass - else: - for month in range(12): - # Mask out the months we don't care about - if month+1 not in months: - model_data[month,:,:] = ma.masked - # Take time average (this will automatically exclude the masked - # months) - model_data = mean(model_data, axis=0) - # Find the maximum along the latitude axis for each longitude point - model_max = amax(model_data, axis=0) - # Find the latitude indices of these maximum values - model_index = argmax(model_data, axis=0) - # Find the latitude values at these indices - model_loc = zeros(size(model_index)) - for i in range(size(model_index)): - model_loc[i] = lat[model_index[i]] - # Use a thicker line for the multi-model mean - if model_name == 'MMM': - lw = 3 - else: - lw = 2 - # Add to plots - ax1.plot(lon, model_max, label=model_name, color=cmip5_colours[j], linewidth=lw) - ax2.plot(lon, model_loc, label=model_name, color=cmip5_colours[j], linewidth=lw) - - # Configure first plot - ax1.set_title('Maximum zonal wind') - ax1.set_xlabel('Longitude') - ax1.set_ylabel('m/s') - ax1.grid(True) - ax1.set_xlim([0, 360]) - # Move the plot over so the legend fits - box = ax1.get_position() - ax1.set_position([box.x0, box.y0, box.width*0.8, box.height]) - fontP = FontProperties() - fontP.set_size('small') - ax1.legend(loc='center left', bbox_to_anchor=(1,0.5), prop=fontP) - - # Configure second plot - ax2.set_title('Latitude of maximum zonal wind') - ax2.set_xlabel('Longitude') - ax2.set_ylabel('Latitude') - ax2.grid(True) - ax2.set_xlim([0, 360]) - box = ax2.get_position() - ax2.set_position([box.x0, box.y0, box.width*0.8, box.height]) - ax2.legend(loc='center left', bbox_to_anchor=(1,0.5), prop=fontP) - - if save: - fig1.savefig(fig_names[0]) - fig2.savefig(fig_names[1]) - else: - fig1.show() - fig2.show() - - -# Command-line interface -if __name__ == "__main__": - - # Make a list of all valid model names - all_models = build_model_list() - all_model_names = [] - for model in all_models: - all_model_names.append(model.name) - all_model_names.append('MMM') - model_names = [] - - model_name = raw_input("Model name: ") - # Check for validity of model name - if model_name in all_model_names: - model_names.append(model_name) - else: - print 'No such model' - while True: - # Keep asking for more models until the user is finished - action = raw_input("Add another model (y/n)? ") - if action == 'y': - model_name = raw_input("Model name: ") - if model_name in all_model_names: - model_names.append(model_name) - else: - print 'No such model' - else: - break - - season = raw_input("Season (djf, mam, jja, son, or annual): ") - - action = raw_input("Save figures (s) or display in window (d)? ") - if action == 's': - save = True - # Get file names for both figures - fig_name1 = raw_input("File name for first figure (maximum Uwind): ") - fig_name2 = raw_input("File name for second figure (latitude of maximum Uwind): ") - fig_names = [fig_name1, fig_name2] - elif action == 'd': - save = False - fig_names = None - - # Make the plots - cmip5_max_uwind(model_names, season, save, fig_names) - - # Repeat until the user wants to exit - while True: - repeat = raw_input("Make another plot (y/n)? ") - if repeat == 'y': - while True: - # Keep asking for changes to parameters until the user is - # finished - changes = raw_input("Enter a paramter to change: (1) models, (2) season, (3) save/display; or enter to continue: ") - if len(changes) == 0: - break - else: - if int(changes) == 1: - # New list of models - model_names = [] - model_name = raw_input("Model name: ") - if model_name in all_model_names: - model_names.append(model_name) - else: - print 'No such model' - while True: - action = raw_input("Add another model (y/n)? ") - if action == 'y': - model_name = raw_input("Model name: ") - if model_name in all_model_names: - model_names.append(model_name) - else: - print 'No such model' - else: - break - elif int(changes) == 2: - # New season - season = raw_input("Season (djf, mam, jja, son, or annual): ") - elif int(changes) == 3: - # Switch from save to display, or vice versa - save = not save - if save: - # Get new figure names - fig_name1 = raw_input("File name for first figure (maximum Uwind): ") - fig_name2 = raw_input("File name for second figure (latitude of maximum Uwind): ") - fig_names = [fig_name1, fig_name2] - # Make the plot - cmip5_max_uwind(model_names, season, save, fig_names) - else: - break - - - - - - diff --git a/cmip5_ocean_climatology_netcdf.py b/cmip5_ocean_climatology_netcdf.py deleted file mode 100644 index 61de093..0000000 --- a/cmip5_ocean_climatology_netcdf.py +++ /dev/null @@ -1,370 +0,0 @@ -from cmip5_paths import * -from cmip5_field import * -from ecco2_field import * -from numpy import * -from netCDF4 import Dataset -from scipy.interpolate import RegularGridInterpolator, griddata -from scipy.spatial import KDTree - -# NB for raijin users: RegularGridInterpolator needs python/2.7.6 but the -# default is 2.7.3. Before running this script, switch them as follows: -# module unload python/2.7.3 -# module unload python/2.7.3-matplotlib -# module load python/2.7.6 -# module load python/2.7.6-matplotlib - -# For the given CMIP5 model, calculates the monthly climatology from 1992-2005 -# inclusive for 4 ocean variables. Interpolates to the ECCO2 grid at the -# northern boundary of ROMS (currently 30S, see ecco2_climatology_netcdf.py) -# and saves to a NetCDF file. Note that to run this script, you must have -# previously run ecco2_climatology_netcdf.py. -# Input: -# model_name = name of model (this must match the list in cmip5_paths.py) -def cmip5_ocean_climatology_netcdf (model_name): - - # Experiment name - expt = 'historical' - # Years over which to calculate climatology - start_year = 1992 - end_year = 2005 - # Northern boundary of ROMS - nbdry = -30.0 - # CMIP5 variable names - var_names = ['thetao', 'so', 'uo', 'vo'] - # Variable names to use in NetCDF file - var_names_output = ['temp', 'salt', 'u', 'v'] - # Units of final variables (note some conversions in cmip5_field) - var_units = ['degC', 'psu', 'm/s', 'm/s'] - # Path to output NetCDF file - output_file = '/short/y99/kaa561/CMIP5_forcing/ocean/' + model_name + '.nc' - # Path to corresponding ECCO2 file (created using - # ecco2_climatology_netcdf.py) - ecco2_file = '/short/y99/kaa561/CMIP5_forcing/ocean/ECCO2.nc' - - # Read ECCO2 grid - id = Dataset(ecco2_file, 'r') - ecco_lon = id.variables['longitude'][:] - ecco_depth = id.variables['depth'][:] - # Save the land mask - mask = id.variables['temp'][0,:,:].mask - id.close() - - # Build a list of Model objects for CMIP5 - all_models = build_model_list() - # Build a corresponding list of all CMIP5 model names - all_model_names = [] - for model in all_models: - all_model_names.append(model.name) - # Find index of model_name in all_model_names, and select the Model object - # at the same index of all_models - # Now model_name has a corresponding Model object - model = all_models[all_model_names.index(model_name)] - - # Loop over variables - for i in range(len(var_names)): - var = var_names[i] - print 'Processing variable ' + var - - # Read monthly climatology for this variable - model_data, model_lon, model_lat, model_depth = cmip5_field(model, expt, var, start_year, end_year) - # Set up array for climatology interpolated to ECCO2 grid - model_data_interp = ma.empty([12, size(ecco_depth), size(ecco_lon)]) - if model_data is not None: - # Interpolate one month at a time - for t in range(size(model_data,0)): - tmp = interp_model2ecco(model_data[t,:,:,:], model_lon, model_lat, model_depth, ecco_lon, nbdry, ecco_depth) - # Mask with ECCO2 land mask - model_data_interp[t,:,:] = ma.masked_where(mask, tmp) - else: - # No data (missing variable in CMIP5 archive) - model_data_interp[:,:,:] = ma.masked - if i == 0: - # Set up NetCDF file on the first iteration - print 'Setting up ' + output_file - id = Dataset(output_file, 'w') - # Define dimensions - id.createDimension('longitude', size(ecco_lon)) - id.createDimension('depth', size(ecco_depth)) - id.createDimension('time', 12) - # Define dimension variables and fill with data - id.createVariable('longitude', 'f8', ('longitude')) - id.variables['longitude'].units = 'degrees' - id.variables['longitude'][:] = ecco_lon - id.createVariable('depth', 'f8', ('depth')) - id.variables['depth'].units = 'metres' - id.variables['depth'][:] = ecco_depth - id.createVariable('time', 'i4', ('time')) - id.variables['time'].units = 'month' - id.variables['time'][:] = arange(1, 12+1) - # Define new variable and fill with data - id.createVariable(var_names_output[i], 'f8', ('time', 'depth', 'longitude')) - id.variables[var_names_output[i]].units = var_units[i] - id.variables[var_names_output[i]][:,:,:] = model_data_interp - id.close() - - -# Interpolate a given CMIP5 field to the ECCO2 grid at the northern boundary -# of ROMS. -# model_data = 3D array (size mxnxo) of model data for the given variable -# model_lon = 1D array (length o) or 2D array (size nxo) of longitude for this -# model -# model_lat = 1D array (length n) or 2D array (size nxo) of latitude for this -# model -# model_depth = 1D array (length m) or 3D array (size mxnxo) of depth for this -# model, in z-coordinates -# ecco_lon = 1D array (length q) of longitude on the ECCO2 grid -# nbdry = latitude of the northern boundary of ROMS -# ecco_depth = 1D array (length p) of depth on the ECCO2 grid -# Output: -# data_interp = 2D array (size pxq) of model data interpolated to the ECCO2 -# grid at the northern boundary of ROMS -def interp_model2ecco (model_data, model_lon, model_lat, model_depth, ecco_lon, nbdry, ecco_depth): - - # Make sure the model's longitude goes from 0 to 360, not -180 to 180 - index = model_lon < 0 - model_lon[index] = model_lon[index] + 360 - - # Check for model depth axis which starts at the bottom - if len(model_depth.shape) == 1: - if model_depth[0] > model_depth[1]: - # Flip it, and the data, so it starts at the surface - model_depth = flipud(model_depth) - model_data = flipud(model_data) - - # Interpolate to northern boundary - - if len(model_lon.shape) == 1 and len(model_lat.shape) == 1: - # Case 1: longitude and latitude are 1D axis variables - - # Find jS and jN, the indices immediately south and immediately north - # of nbdry - if model_lat[0] > model_lat[1]: - # Latitude is decreasing - # Find the first index south of nbdry - jS = nonzero(model_lat < nbdry)[0][0] - # The index before it is the last index north of nbdry - jN = jS - 1 - elif model_lat[0] < model_lat[1]: - # Latitude is increasing - # Find the first index north of nbdry - jN = nonzero(model_lat > nbdry)[0][0] - # The index before it is the last index south of nbdry - jS = jN - 1 - # Linearly interpolate model_data to nbdry - model_dataS = model_data[:,jS,:] - model_dataN = model_data[:,jN,:] - model_data_bdry = (model_dataN - model_dataS)/(model_lat[jN] - model_lat[jS])*(nbdry - model_lat[jS]) + model_dataS - - elif len(model_lon.shape) == 2 and len(model_lat.shape) == 2: - # Case 2: longitude and latitude are 2D variables, i.e. rotated grid - - # Set up empty arrays for data and longitude interpolated to nbdry - model_data_bdry = ma.empty((size(model_data,0), size(model_data,2))) - model_lon_bdry = ma.empty((size(model_data,2))) - if len(model_depth.shape) == 3: - # We will also have to interpolate depth to nbdry - model_depth_bdry = ma.empty(shape(model_data_bdry)) - # Loop over longitude axis and interpolate one index at a time - for i in range(size(model_data,2)): - # Select current latitude values - model_lat_tmp = model_lat[:,i] - # Check for non-monotonic latitude values - dlat = model_lat_tmp[1:] - model_lat_tmp[0:-1] - if any(dlat < 0) and any(dlat > 0): - # Latitude doesn't strictly increase or strictly decrease with j - # A few models have rotated grids where latitude is monotonic - # in the middle but changes in the other direction near the ends - # We will just mask out the ends, one index at a time, until - # the leftover values in the middle are monotonic - # This works because the non-monotonic sections are near the - # poles, and far away from nbdry (currently 30S) - model_lat_tmp = ma.masked_array(model_lat_tmp) - posn = 1 - while True: - model_lat_tmp[posn-1] = ma.masked - model_lat_tmp[-posn] = ma.masked - dlat = model_lat_tmp[1:] - model_lat_tmp[0:-1] - if all(dlat < 0) or all(dlat > 0): - break - posn += 1 - else: - posn = 0 - if model_lat_tmp[posn] > model_lat_tmp[posn+1]: - # Latitude is decreasing - # Find the first index south of nbdry - jS = nonzero(model_lat_tmp < nbdry)[0][0] - # The index before it is the last index north of nbdry - jN = jS - 1 - elif model_lat_tmp[posn] < model_lat_tmp[posn+1]: - # Latitude is increasing - # Find the first index north of nbdry - jN = nonzero(model_lat_tmp > nbdry)[0][0] - # The index before it is the last index south of nbdry - jS = jN - 1 - # Linearly interpolate model_data to nbdry - model_dataS = model_data[:,jS,i] - model_dataN = model_data[:,jN,i] - model_data_bdry[:,i] = (model_dataN - model_dataS)/(model_lat[jN,i] - model_lat[jS,i])*(nbdry - model_lat[jS,i]) + model_dataS - # Linearly interpolate model_lon to nbdry - model_lonS = model_lon[jS,i] - model_lonN = model_lon[jN,i] - if model_lonS < 1 and model_lonN > 300: - # The jump in longitude from almost-360 to almost-0 happens - # between lonS and lonN. Take mod 360 to fix this. - model_lonS += 360 - model_lon_bdry[i] = (model_lonN - model_lonS)/(model_lat[jN,i] - model_lat[jS,i])*(nbdry - model_lat[jS,i]) + model_lonS - if len(model_depth.shape) == 3: - # Linearly interpolate depth to nbdry - model_depthS = model_depth[:,jS,i] - model_depthN = model_depth[:,jN,i] - model_depth_bdry[:,i] = (model_depthN - model_depthS)/(model_lat[jN,i] - model_lat[jS,i])*(nbdry - model_lat[jS,i]) + model_depthS - # Overwrite model_lon with the interpolated values; now it's 1D - model_lon = model_lon_bdry - if len(model_depth.shape) == 3: - # Overwrite model_depth with the interpolated values; now it's 2D - model_depth = model_depth_bdry - - if model_lon[-2] == model_lon[0] and model_lon[-1] == model_lon[1]: - # A few models have grids which overlap by 2 indices at the periodic - # boundary; delete these last 2 indices - model_lon = model_lon[0:-2] - model_data_bdry = model_data_bdry[:,0:-2] - - # Most grids will have longitude jump from almost-360 to almost-0 - # at some point in the interior. Find the index where this happens, and - # rearrange the model data and axes so the jump is moved to the periodic - # boundary. This way model_lon will be monotonic, which is needed for - # interpolation. - dlon = model_lon[1:] - model_lon[0:-1] - index = argmax(abs(dlon)) - if dlon[index] < -300: - # Add 1 to index because dlon is one index shorter than model_lon - jump = index+1 - # Rearrange longitude values and data - model_lon_new = zeros(size(model_lon)) - model_lon_new[0:-jump] = model_lon[jump:] - model_lon_new[-jump:] = model_lon[0:jump] - model_lon = model_lon_new - model_data_new = ma.empty(shape(model_data_bdry)) - model_data_new[:,0:-jump] = model_data_bdry[:,jump:] - model_data_new[:,-jump:] = model_data_bdry[:,0:jump] - model_data_bdry = model_data_new - if len(model_depth.shape) == 2: - # Rearrange depth values - model_depth_new = ma.empty(shape(model_depth)) - model_depth_new[:,0:-jump] = model_depth[:,jump:] - model_depth_new[:,-jump:] = model_depth[:,0:jump] - model_depth = model_depth_new - - # CMIP5 model axes don't wrap around; there is a gap between almost-180W - # and almost-180E (these are 0 to 360 but you get the idea) and depending - # on the grid, we may need to interpolate in this gap. - # So copy the last longitude value (mod 360) to the beginning, and the - # first longitude value (mod 360) to the end. - model_lon_wrap = zeros(size(model_lon)+2) - model_lon_wrap[0] = model_lon[-1] - 360 - model_lon_wrap[1:-1] = model_lon - model_lon_wrap[-1] = model_lon[0] + 360 - model_lon = model_lon_wrap - # Copy the westernmost and easternmost data points to match - model_data_wrap = ma.empty((size(model_data_bdry,0), size(model_lon))) - model_data_wrap[:,1:-1] = model_data_bdry - model_data_wrap[:,0] = model_data_bdry[:,-1] - model_data_wrap[:,-1] = model_data_bdry[:,0] - model_data_bdry = model_data_wrap - if len(model_depth.shape) == 2: - # Copy the westernmost and easternmost depth values to match - model_depth_wrap = ma.empty((size(model_depth,0), size(model_lon))) - model_depth_wrap[:,1:-1] = model_depth - model_depth_wrap[:,0] = model_depth[:,-1] - model_depth_wrap[:,-1] = model_depth[:,0] - model_depth = model_depth_wrap - - if amin(model_depth) >= amin(ecco_depth): - # The model grid doesn't extend far enough to the surface to be able - # to interpolate to ECCO2's surface points (5 metres depth). - # Add another point at 0m, and just copy the data values closest - # to the surface. - if len(model_depth.shape) == 1: - model_depth_new = zeros(size(model_depth)+1) - model_depth_new[1:] = model_depth - elif len(model_depth.shape) == 2: - model_depth_new = zeros(size(model_depth,0)+1, size(model_depth,1)) - model_depth_new[1:,:] = model_depth - model_depth = model_depth_new - model_data_new = ma.empty((size(model_data_bdry,0)+1, size(model_lon))) - model_data_new[1:,:] = model_data_bdry - model_data_new[0,:] = model_data_new[1,:] - model_data_bdry = model_data_new - if amax(model_depth) <= amax(ecco_depth): - # The model grid doesn't extend deep enough to be able to interpolate - # to ECCO2's bottom points (almost 6000m depth). - # Add another point at 6000m, and just copy the data values closest - # to the bottom. - if len(model_depth.shape) == 1: - model_depth_new = zeros(size(model_depth)+1) - model_depth_new[0:-1] = model_depth - model_depth_new[-1] = 6000.0 - elif len(model_depth.shape) == 2: - model_depth_new = zeros((size(model_depth,0)+1, size(model_depth,1))) - model_depth_new[0:-1,:] = model_depth - model_depth_new[-1,:] = 6000.0 - model_depth = model_depth_new - model_data_new = ma.empty((size(model_data_bdry,0)+1, size(model_lon))) - model_data_new[0:-1,:] = model_data_bdry - model_data_new[-1,:] = model_data_new[-2,:] - model_data_bdry = model_data_new - - # Fill in masked values with nearest neighbours so they don't screw up - # the interpolation - # I got this code from Stack Exchange, not really sure how it works - j,i = mgrid[0:model_data_bdry.shape[0], 0:model_data_bdry.shape[1]] - jigood = array((j[~model_data_bdry.mask], i[~model_data_bdry.mask])).T - jibad = array((j[model_data_bdry.mask], i[model_data_bdry.mask])).T - model_data_bdry[model_data_bdry.mask] = model_data_bdry[~model_data_bdry.mask][KDTree(jigood).query(jibad)[1]] - - # Get a 2D mesh of ECCO2 longitude and depth - ecco_lon_2d, ecco_depth_2d = meshgrid(ecco_lon, ecco_depth) - - if len(model_depth.shape) == 1: - - # Build an interpolation function for model_data - interp_function = RegularGridInterpolator((model_depth, model_lon), model_data_bdry) - # Call it for the ECCO2 grid - data_interp = interp_function((ecco_depth_2d, ecco_lon_2d)) - - elif len(model_depth.shape) == 2: - - # Since the grid isn't regular in lon-depth space, we have to use - # griddata instead of RegularGridInterpolator - # First copy the model's longitude values into a 2D array of dimension - # depth x longitude - model_lon_2d = tile(model_lon, (size(model_data_bdry,0), 1)) - # Now set up an nx2 array containing the coordinates of each point - # in the flattened array of model data - points = empty([size(model_data_bdry),2]) - points[:,0] = ravel(model_lon_2d) - points[:,1] = ravel(model_depth) - # Also flatten the model data - values = ravel(model_data_bdry) - # Now set up an mx2 array containing the coordinates of each point - # we want to interpolate to, in the ECCO2 grid - xi = empty([size(ecco_lon_2d),2]) - xi[:,0] = ravel(ecco_lon_2d) - xi[:,1] = ravel(ecco_depth_2d) - # Now call griddata - result = griddata(points, values, xi) - # Un-flatten the result - data_interp = reshape(result, shape(ecco_lon_2d)) - - return data_interp - - -# Command-line interface -if __name__ == "__main__": - - # Process one model at a time - models = build_model_list() - for model in models: - print model.name - cmip5_ocean_climatology_netcdf(model.name) diff --git a/cmip5_paths.py b/cmip5_paths.py deleted file mode 100644 index 985306e..0000000 --- a/cmip5_paths.py +++ /dev/null @@ -1,94 +0,0 @@ -from os.path import exists - -# Model object containing name of CMIP5 model, institution and ensemble member. -class Model: - - def __init__ (self, inst, name): - - self.name = name - self.inst = inst - self.ens = 'r1i1p1' - if self.name == 'CESM1-WACCM': - self.ens = 'r2i1p1' - - - # Return the directory containing monthly output for the given variable and - # experiment. - # Input: - # expt = string containing name of experiment, eg 'historical' - # var_name = string containing name of variable, eg 'tas' - # Output: dir = string containing directory where the model output is - # stored. If this model output doesn't exist, return an - # empty string. - def get_directory (self, expt, var_name): - - # Figure out whether it is an atmosphere or ocean variable - if var_name in ['ps', 'tas', 'huss', 'clt', 'uas', 'vas', 'pr', 'prsn', 'evspsbl', 'rsds', 'rlds']: - realm = 'atmos' - elif var_name in ['thetao', 'so', 'uo', 'vo', 'zos']: - realm = 'ocean' - else: - print 'Unknown variable' - # Exit early - return '' - - # Build typical directory structure in ua6's archive on raijin - dir = '/g/data1/ua6/drstree/CMIP5/GCM/' + self.inst + '/' + self.name + '/' + expt + '/mon/' + realm + '/' + var_name + '/' + self.ens + '/' - - # Exceptions - if self.name == 'GFDL-ESM2M' and expt == 'rcp85' and var_name == 'clt': - dir = '/g/data/ua6/unofficial-ESG-replica/tmp/tree/esgdata.gfdl.noaa.gov/thredds/fileServer/gfdl_dataroot/NOAA-GFDL/GFDL-ESM2M/rcp85/mon/atmos/Amon/r1i1p1/v20110601/clt/' - - # Check if this directory actually exists, or if the data is missing. - if not exists(dir): - # Return an empty string - dir = '' - - return dir - - -# Build an array of Model objects for the 39 CMIP5 models used in this project. -def build_model_list (): - - Models = [] - Models.append(Model('BCC', 'bcc-csm1-1')) - Models.append(Model('BCC', 'bcc-csm1-1-m')) - Models.append(Model('BNU', 'BNU-ESM')) - Models.append(Model('CCCMA', 'CanESM2')) - Models.append(Model('CMCC', 'CMCC-CM')) - Models.append(Model('CMCC', 'CMCC-CMS')) - Models.append(Model('CNRM', 'CNRM-CM5')) - Models.append(Model('CSIRO-BOM', 'ACCESS1-0')) - Models.append(Model('CSIRO-BOM', 'ACCESS1-3')) - Models.append(Model('CSIRO-QCCCE', 'CSIRO-Mk3-6-0')) - Models.append(Model('FIO', 'FIO-ESM')) - Models.append(Model('ICHEC', 'EC-EARTH')) - Models.append(Model('INM', 'inmcm4')) - Models.append(Model('IPSL', 'IPSL-CM5A-LR')) - Models.append(Model('IPSL', 'IPSL-CM5A-MR')) - Models.append(Model('IPSL', 'IPSL-CM5B-LR')) - Models.append(Model('LASG-CESS', 'FGOALS-g2')) - Models.append(Model('MIROC', 'MIROC-ESM')) - Models.append(Model('MIROC', 'MIROC-ESM-CHEM')) - Models.append(Model('MIROC', 'MIROC5')) - Models.append(Model('MOHC', 'HadGEM2-CC')) - Models.append(Model('MOHC', 'HadGEM2-ES')) - Models.append(Model('MPI-M', 'MPI-ESM-LR')) - Models.append(Model('MPI-M', 'MPI-ESM-MR')) - Models.append(Model('MRI', 'MRI-CGCM3')) - Models.append(Model('NASA-GISS', 'GISS-E2-H')) - Models.append(Model('NASA-GISS', 'GISS-E2-H-CC')) - Models.append(Model('NASA-GISS', 'GISS-E2-R')) - Models.append(Model('NASA-GISS', 'GISS-E2-R-CC')) - Models.append(Model('NCAR', 'CCSM4')) - Models.append(Model('NCC', 'NorESM1-M')) - Models.append(Model('NCC', 'NorESM1-ME')) - Models.append(Model('NIMR-KMA', 'HadGEM2-AO')) - Models.append(Model('NOAA-GFDL', 'GFDL-CM3')) - Models.append(Model('NOAA-GFDL', 'GFDL-ESM2G')) - Models.append(Model('NOAA-GFDL', 'GFDL-ESM2M')) - Models.append(Model('NSF-DOE-NCAR', 'CESM1-BGC')) - Models.append(Model('NSF-DOE-NCAR', 'CESM1-CAM5')) - Models.append(Model('NSF-DOE-NCAR', 'CESM1-WACCM')) - - return Models diff --git a/ecco2_climatology_netcdf.py b/ecco2_climatology_netcdf.py deleted file mode 100644 index 58471cd..0000000 --- a/ecco2_climatology_netcdf.py +++ /dev/null @@ -1,67 +0,0 @@ -from ecco2_field import * -from numpy import * -from netCDF4 import Dataset - -# Calculate the ECCO2 monthly climatology from 1992-2005 inclusive, interpolated -# to the northern boundary of the ROMS domain (currently 30S, this is set in -# ecco2_field.py), for 4 ocean variables. Save to a NetCDF file. -def ecco2_climatology_netcdf (): - - # Date range for climatology - start_year = 1992 - end_year = 2005 - # Names of ECCo2 variables - var_names = ['THETA', 'SALT', 'UVEL', 'VVEL'] - # Variable names to use for NetCDF file - var_names_output = ['temp', 'salt', 'u', 'v'] - # Units of final variables (note there are some conversions in ecco2_field) - var_units = ['degC', 'psu', 'm/s', 'm/s'] - # Path to output NetCDF file - output_file = '/short/y99/kaa561/CMIP5_forcing/ocean/ECCO2.nc' - - # Loop over variables - for i in range(len(var_names)): - var = var_names[i] - print 'Processing variable ' + var - # Read data and grid - data, lon, depth = ecco2_field(var, start_year, end_year) - if i == 0: - # Save mask of the first variable (temperature) - orig_mask = data.mask - # Create NetCDF file on the first iteration - print 'Setting up ' + output_file - id = Dataset(output_file, 'w') - # Define dimensions - id.createDimension('longitude', size(lon)) - id.createDimension('depth', size(depth)) - id.createDimension('time', 12) - # Define dimension variables and fill with data - id.createVariable('longitude', 'f8', ('longitude')) - id.variables['longitude'].units = 'degrees' - id.variables['longitude'][:] = lon - id.createVariable('depth', 'f8', ('depth')) - id.variables['depth'].units = 'metres' - id.variables['depth'][:] = depth - id.createVariable('time', 'i4', ('time')) - id.variables['time'].units = 'month' - id.variables['time'][:] = arange(1, 12+1) - if all(~data.mask): - # There is no mask on this data - # This happens for u and v which are filled with zeros in the mask - # Apply the mask from temperature, saved earlier - data = ma.masked_where(orig_mask, data) - # Define a new variable in the NetCDF file and fill with data - id.createVariable(var_names_output[i], 'f8', ('time', 'depth', 'longitude')) - id.variables[var_names_output[i]].units = var_units[i] - id.variables[var_names_output[i]][:,:,:] = data - id.close() - - -# Command-line interface -if __name__ == "__main__": - - ecco2_climatology_netcdf() - - - - diff --git a/ecco2_field.py b/ecco2_field.py deleted file mode 100644 index 73dd2a3..0000000 --- a/ecco2_field.py +++ /dev/null @@ -1,71 +0,0 @@ -from numpy import * -from netCDF4 import Dataset - -# Read ECCO2 data for the given variable name, between the given start and end -# years. Interpolate to the northern boundary of the circumpolar ROMS domain. -# Return the monthly climatology. -# Input: -# var_name = string containing name of variable, eg 'THETA' -# start_year, end_year = integers containing years over which to calculate the -# monthly climatology, from the beginning of start_year -# to the end of end_year. Therefore, if start_year = -# end year, this script will read one year of output -# with no climatological averaging. -# Output: -# ecco2_data = 3D array of ECCO2 data, with dimension month x depth x longitude -# ecco2_lon = 1D array containing longitude values -# ecco2_depth = 1D array containing depth values -def ecco2_field (var_name, start_year, end_year): - - # Latitude of the northern boundary of the circumpolar ROMS domain - nbdry = -30 - # Directory where ECCO2 data is stored - ecco2_dir = '/short/m68/kaa561/ROMS-CICE-MCT/data/ECCO2/raw/' - # Middle of ECCO2 file names - ecco2_mid = '.1440x720x50.' - - # Read ECCO2 latitude, longitude, and depth from the first file - id = Dataset(ecco2_dir + 'THETA' + ecco2_mid + str(start_year) + '01.nc', 'r') - ecco2_lat = id.variables['LATITUDE_T'][:] - ecco2_lon = id.variables['LONGITUDE_T'][:] - ecco2_depth = id.variables['DEPTH_T'][:] - id.close() - # Find the first index north of nbdry, and subtract 1 to find the last - # index south of nbdry - j_max = nonzero(ecco2_lat > nbdry)[0][0] - j_min = j_max - 1 - # Only save the ECCO2 latitude values at these indices - ecco2_lat = ecco2_lat[j_min:j_max+1] - - ecco2_data = None - # Loop over years and months - for year in range(start_year, end_year+1): - for month in range(12): - # Construct filename - if (month+1) < 10: - month_str = '0' + str(month+1) - else: - month_str = str(month+1) - ecco2_file = ecco2_dir + str(var_name) + ecco2_mid + str(year) + month_str + '.nc' - # Read data - id = Dataset(ecco2_file, 'r') - data = id.variables[var_name][0,:,j_min:j_max+1,:] - id.close() - - # Linearly interpolate to nbdry - data_interp = (data[:,1,:]-data[:,0,:])/(ecco2_lat[1]-ecco2_lat[0])*(nbdry-ecco2_lat[0]) + data[:,0,:] - if ecco2_data is None: - # Create empty array of dimension time x depth x longitude - ecco2_data = ma.empty([12, size(ecco2_depth), size(ecco2_lon)]) - # Initialise with zeros and mask with land mask - for t in range(12): - ecco2_data[t,:,:] = data_interp[:,:]*0.0 - # Add to master array - ecco2_data[month,:,:] += data_interp[:,:] - - # Divide master array by number of years to convert from monthly sums to - # monthly averages - ecco2_data /= (end_year-start_year+1) - - return ecco2_data, ecco2_lon, ecco2_depth - diff --git a/eraint_climatology_netcdf.py b/eraint_climatology_netcdf.py deleted file mode 100644 index cf24366..0000000 --- a/eraint_climatology_netcdf.py +++ /dev/null @@ -1,56 +0,0 @@ -from eraint_field import * -from numpy import * -from netCDF4 import Dataset - -# Calculate the ERA-Interim monthly climatology from 1992-2005 inclusive, for 11 -# atmospheric variables (all the variables which ROMS and/or FESOM depend on). -# Save to a NetCDF file. -def eraint_climatology_netcdf (): - - # Date range for climatology - start_year = 1992 - end_year = 2005 - # Names of ERA-Interim variables - var_names = ['sp', 't2m', 'd2m', 'tcc', 'u10', 'v10', 'tp', 'sf', 'e', 'ssrd', 'strd'] - # Variable names to use for NetCDF file - var_names_output = ['Pair', 'Tair', 'Hair', 'cloud', 'Uwind', 'Vwind', 'precip', 'snow', 'evap', 'swrad', 'lwrad'] - # Units of final variables (note there are some conversions in eraint_field) - var_units = ['kPa', 'degC', '1', '%', 'm/s', 'm/s', '10^6 kg/m^2/s', '10^6 kg/m^2/s', '10^6 kg/m^2/s', 'W/m^2', 'W/m^2'] - # Path to output NetCDF file - output_file = '/short/y99/kaa561/CMIP5_forcing/atmos/ERA-Interim.nc' - - # Loop over variables - for i in range(len(var_names)): - var = var_names[i] - print 'Processing variable ' + var - # Read data and grid - data, lon, lat = eraint_field(var, start_year, end_year) - if i == 0: - # Create NetCDF file on the first iteration - print 'Setting up ' + output_file - id = Dataset(output_file, 'w') - # Define dimensions - id.createDimension('longitude', size(lon)) - id.createDimension('latitude', size(lat)) - id.createDimension('time', 12) - # Define dimension variables and fill with data - 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('time', 'f8', ('time')) - id.variables['time'].units = 'month' - id.variables['time'][:] = arange(1, 12+1) - # Define a new variable in the NetCDF file and fill with data - id.createVariable(var_names_output[i], 'f8', ('time', 'latitude', 'longitude')) - id.variables[var_names_output[i]].units = var_units[i] - id.variables[var_names_output[i]][:,:,:] = data - id.close() - - -# Command-line interface -if __name__ == "__main__": - - eraint_climatology_netcdf() diff --git a/eraint_field.py b/eraint_field.py deleted file mode 100644 index 76148d9..0000000 --- a/eraint_field.py +++ /dev/null @@ -1,99 +0,0 @@ -from numpy import * -from netCDF4 import Dataset - -# Read ERA-Interim monthly data for the given variable name, between the given -# start and end years. Return the monthly climatology. -# Input: -# var_name = string containing name of variable, eg 't2m' -# start_year, end_year = integers containing years over which to calculate the -# monthly climatology, from the beginning of start_year -# to the end of end_year. Therefore, if start_year = -# end year, this script will read one year of output -# with no climatological averaging. -# Output: -# era_data = 3D array of ERA-Interim data, with dimension month x latitude x -# longitude, possibly with units converted to be more comparable -# to CMIP5 models -# era_lon = 1D array containing longitude values -# era_lat = 1D array containing latitude values -def eraint_field (var_name, start_year, end_year): - - # Directory where ERA-Interim monthly averaged data is stored - era_dir = '/short/y99/kaa561/FESOM/ERA_Interim_monthly/' - # String that ERA-Interim files end with - era_tail = '_monthly_orig.nc' - - # Latent heat of vapourisation, J/kg - Lv = 2.5e6 - # Ideal gas constant for water vapour, J/K/kg - Rv = 461.5 - # Density of water, kg/m^3 - rho_w = 1e3 - # Conversion from K to C - degKtoC = -273.15 - - # Read ERA-Interim latitude and longitude - id = Dataset(era_dir + 'AN_' + str(start_year) + era_tail, 'r') - era_lat = id.variables['lat'][:] - era_lon = id.variables['lon'][:] - id.close() - - # Figure out how ERA-Interim filename will start; the atmospheric - # variables for each year are split between 3 different files - if var_name in ['sp', 't2m', 'd2m', 'tcc', 'u10', 'v10']: - era_head = era_dir + 'AN_' - elif var_name in ['tp', 'sf']: - era_head = era_dir + 'FC_' - elif var_name in ['e', 'ssrd', 'strd']: - era_head = era_dir + 'ER_' - - era_data = None - # Loop over years - for year in range(start_year, end_year+1): - - # Construct filename - era_file = era_head + str(year) + era_tail - # Read data - id = Dataset(era_file, 'r') - data = id.variables[var_name][:,:,:] - - # Perform conversions if necessary - if var_name == 'sp': - # Convert pressure from Pa to kPa - data = 1e-3*data - elif var_name == 't2m': - # Convert temperature from K to C - data = data + degKtoC - elif var_name == 'd2m': - # Calculate specific humidity from dewpoint temperature - # and surface pressure - d2m = data - sp = id.variables['sp'][:,:,:] - # Intermediate step to calculate vapour pressure - e = 611*exp(Lv/Rv*(1/273.0 - 1/d2m)) - data = 0.622*e/(sp - 0.378*e) - elif var_name == 'tcc': - # Convert total cloud cover from fraction to percent - data = data*100 - elif var_name in ['tp', 'sf', 'e']: - # Convert precip/snowfall/evap from kg/m^2/s to 1e-6 kg/m^2/s - data = 1e6*data*rho_w/(12*60*60) - if var_name == 'e': - # ERA-Interim defines evaporation as negative; fix this - data = -1*data - elif var_name in ['ssrd', 'strd']: - # Convert from J/m^2, integrated over 12 hours, to W/m^2 - data = data/(12*60*60) - id.close() - - # Add to master array - if era_data is None: - era_data = data - else: - era_data += data - - # Divide master array by number of years to convert from monthly sums to - # monthly averages - era_data /= (end_year-start_year+1) - - return era_data, era_lon, era_lat diff --git a/file_guide.txt b/file_guide.txt index 2a4b472..a86acc9 100644 --- a/file_guide.txt +++ b/file_guide.txt @@ -326,155 +326,6 @@ timeseries_3D.py: Calculates and plots timeseries of total ocean heat content, or averages file, and the log file. -***CMIP5 ANALYSIS*** - -cmip5_paths.py: Builds an array of Model objects, one for each of the 39 CMIP5 - models used in Kaitlin's project. Model objects contain - information about the model, and most usefully a class function - get_directory which returns the directory containing monthly - output for the specified variable and experiment. - To run: The function build_model_list is designed to be called - by other functions; see for example - cmip5_eraint_rms_errors.py. Note that for the class - function get_directory to work, you must be running on - raijin and be part of the ua6 project. - -eraint_field.py: Reads ERA-Interim monthly data for the given variable name, - between the given start and end years. Returns the monthly - climatology. - To run: The function eraint_field is designed to be called - by other functions; see for example - eraint_climatology_netcdf.py. If you are using it, - be sure to edit the parameters near the top of the - file (eg paths to ERA-Interim monthly data). - -ecco2_field.py: Reads ECCO2 data for the given variable name, between the given - start and end years. Interpolates to the northern boundary of - the circumpolar ROMS domain (currently 30S). Returns the monthly - climatology. - To run: The function ecco2_field is designed to be called by - other functions; see for example - ecco2_climatology_netcdf.py. If you are using it, be - sure to edit the parameters near the top of the file - (eg paths to ECCO2 data, and the location of the ROMS - northern boundary). - -cmip5_field.py: Reads CMIP5 output for the given model, experiment, and variable - name. Returns the monthly climatology as well as the grid - (longitude, latitude, and depth arrays). - To run: The function cmip5_field is designed to be called by - other functions; see for example - cmip5_atmos_climatology_netcdf.py. Note that for this - function to work, you must be running on raijin and be - part of the ua6 project. - -eraint_climatology_netcdf.py: Calculates the ERA-Interim monthly climatology - from 1992-2005 inclusive, for 11 atmospheric - variables (all the variables which ROMS and/or - FESOM depend on). Save to a NetCDF file. - To run: First edit the variable output_file, near - the top of the file, to suit your - filesystem. Then open python or ipython - and type - "run eraint_climatology_netcdf.py". - -ecco2_climatology_netcdf.py: Like eraint_climatology_netcdf.py, but for ECCO2 - ocean variables rather than ERA-Interim atmoshperic - variables, interpolated to the northern boundary - of the ROMS domain (set in ecco2_field.py). - To run: First edit the variable output_file, near - the top of the file, to suit your - filesystem. Then open python or ipython and - type "run ecco2_climatology_netcdf.py". - -cmip5_atmos_climatology_netcdf.py: For the given CMIP5 model, calculates the - monthly climatology from 1992-2005 inclusive for 11 - atmospheric variables. Interpolates to the ERA-Interim - grid and saves to a NetCDF file. Note that to run this - script, you must have previously run - eraint_climatology_netcdf.py. - To run: First edit the variables output_file and - eraint_file, near the top of the file, to suit - your filesystem. 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). To call this function - for all CMIP5 models, open python or ipython and - type "run cmip5_atmos_climatology_netcdf.py". - To call this function for just one model, e.g. - ACCESS1-0, open python or ipython and type - "from cmip5_atmos_climatology_netcdf import *" - and then - "cmip5_atmos_climatology_netcdf('ACCESS1-0')". - -cmip5_ocean_climatology_netcdf.py: For the given CMIP5 model, calculates the - monthly climatology from 1992-2005 inclusive for 4 ocean - variables. Interpolates to the ECCO2 grid at the northern - boundary of ROMS and saves to a NetCDF file. Note that to - run this script, you must have previously run - ecco2_climatology_netcdf.py. - To run: First edit the variables output_file and - ecco2_file, near the top of the file, to suit - your filesystem. 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). To call this function for - all CMIP5 models, open python or ipython and type - "run cmip5_ocean_climatology_netcdf.py". - To call this function for just one model, e.g. - ACCESS1-0, open python or ipython and type - "from cmip5_ocean_climatology_netcdf import *" - and then - "cmip5_ocean_climatology_netcdf('ACCESS1-0')". - -mmm_atmos_netcdf.py: Calculate the multi-model mean of atmospheric climatology - files created using cmip5_atmos_climatology_netcdf.py. - To run: First edit the variable "directory" near the top - of the file. Then open python or ipython and type - "run mmm_atmos_netcdf.py". - -mmm_ocean_netcdf.py: Calculate the multi-model mean of ocean climatology files - created using cmip5_ocean_climatology_netcdf.py. - To run: First edit the variable "directory" near the top of - the file. Then open python or ipython and type - "run mmm_ocean_netcdf.py". - -cmip5_eraint_rms_errors.py: Calculates root-mean-square errors (as in Gleckler - et al., 2008) for each of 39 CMIP5 models and the - multi-model mean, with respect to 11 ERA-Interim - atmospheric variables. The domain is the Southern - Ocean (all longitudes, and latitudes from the - northern boundary of ROMS to the southernmost ocean - point not in an ice shelf cavity) and the monthly - climatology averaged over 1992-2005 inclusive. - Also calculate the relative errors as in Gleckler - et al. and make a "portrait plot" of coloured tiles - in a model vs. variable matrix. Save both rms errors - and relative errors into text files. Note that to - run this script, you must have previously run - eraint_climatology_netcdf.py, - cmip5_atmos_climatology_netcdf.py, and - mmm_atmos_netcdf.py. - To run: First edit the variable "directory", near - the top of the file, to suit your - filesystem, and "roms_grid" to point to the - roms grid file. Then open python or ipython - and type "run cmip5_eraint_rms_errors.py". - -cmip5_ecco2_rms_errors.py: Same as cmip5_eraint_rms_errors.py, but for ECCO2 - ocean variables rather than ERA-Interim atmospheric - variables, at the northern boundary of ROMS - (currently 30S). Note that to run this script, you - must have previously run - ecco2_climatology_netcdf.py, - cmip5_ocean_climatology_netcdf.py, and - mmm_ocean_netcdf.py. - To run: First edit the variable "directory", near - the top of the file, to suit your filesystem. - Then open python or ipython and type - "run cmip5_ecco2_rms_errors.py". - - ***NICE FIGURES*** zice.py: Saves a contour plot of ice shelf draft, with the land masked in white, @@ -574,47 +425,6 @@ zonal_plot.py: Creates a depth vs latitude plot of the given variable. You can only have to specify which parameters have changed since the last plot. -cmip5_plot.py: Compare output from CMIP5 models to ERA-Interim (for atmosphere - variables) or ECCO2 (for ocean variables) by plotting the given - variable, time-averaged over the given season and zonally - averaged over the Southern Ocean (for atmosphere variables) or - the northern boundary of ROMS (for ocean variables). The plot - will have the given variable on the x-axis and latitude (for - atmosphere variables) or depth (for ocean variables) on the - y-axis. Note that in order to run this script, you must first - run eraint_climatology_netcdf.py, ecco2_climatology_netcdf.py, - cmip5_atmos_climatology_netcdf.py, - cmip5_ocean_climatology_netcdf.py, mmm_atmos_netcdf.py, and - mmm_ocean_netcdf.py to generate the necessary NetCDF files. - To run: First edit the variable "directory", near the top of the - file, to suit your filesystem, and the variable - "roms_grid" to point to your ROMS grid file. (Note that - your choice of roms_grid will affect the latitude bounds - on atmosphere plots but not ocean plots, which read - NetCDF files which have already been interpolated to - the northern boundary of ROMS). Then open python or - ipython and type "run cmip5_plot.py". The script will - prompt you for the variable name, season code, models - to include ('MMM' is an option for multi-model mean), - and whether to save the plot (and if so, what the - filename should be) or display it on the screen. - As with circumpolar_plot.py, circumpolar_cice_plot.py, - and zonal_plot.py, the interface will repeat as many - times as you like, and you only have to specify which - parameters have changed since the last plot. - -cmip5_all_plots.py: Call cmip5_plot.py for all models (including the multi-model - mean), all variables (atmosphere and ocean), and all - seasons. Note that in order to run this script, you must - have previously run eraint_climatology_netcdf.py, - ecco2_climatology_netcdf.py, - cmip5_atmos_climatology_netcdf.py, - cmip5_ocean_climatology_netcdf.py, mmm_atmos_netcdf.py, and - mmm_ocean_netcdf.py. Also, the directory "cmip5/" must - exist. - To run: Open python or ipython and type - "run cmip5_all_plots.py". - overturning_plot.py: Calculate the meridional overturning streamfunction at the last timestep of the given ROMS history/averages file and make a contour plot in z-space. @@ -623,27 +433,6 @@ overturning_plot.py: Calculate the meridional overturning streamfunction at the for the path to the ROMS history/averages file and the filename with which to save the plot. -cmip5_max_uwind.py: Create two plots: (1) the maximum zonal wind speed between - 30S and 65S, and (2) the latitude of that maximum zonal - wind speed, both against longitude. Plot results for - ERA-Interim as well as the given CMIP5 models, averaged - over the given season, for the 1992-2005 climatology. - Note that in order to run this script, you must have - previously run eraint_climatology_netcdf.py, - cmip5_atmos_climatology_netcdf.py, and mmm_atmos_netcdf.py. - To run: First edit the variable "directory", near the top - of the file, to suit your filesystem. Then open - python or ipython and type "run cmip5_max_uwind.py". - The script will prompt you for the models to include - ('MMM' is an option for multi-model-mean), the - season code, and whether to save the plots (and if - so, what the filenames should be) or display them - on the screen. As with circumpolar_plot.py, - circumpolar_cice_plot.py, zonal_plot.py, and - cmip5_plot.py, the interface will repeat as many - times as you like, and you only have to specify - which parameters have changed since the last plot. - uv_vectorplot.py: Make a circumpolar Antarctic plot of speed overlaid with velocity vectors at the given depth (surface, bottom, or vertically aveaged). @@ -788,6 +577,40 @@ make_zonal_slices.sh: Example of a bash script that calls sose_roms_seasonal.py Then just type "./make_zonal_slices.sh" in the shell and it will go. +aice_hi_seasonal.py: Creates a 4x2 plot of seasonally averaged sea ice + concentration (top row) and thickness (bottom row) over + the last year of simulation. + To run: Open python or ipython and type + "run aice_hi_seasonal.py". The script will prompt + you for the CICE output file and whether you want + to save the figure (and if so, what filename) or + display it on the screen. This script does not + repeat. + +freezingpt_slice.py: Plot the difference from the freezing temperature for a + single i-slice (depth versus y) with no time-averaging, + spatial averaging, or interpolation. This helps to show + where there is spurious supercooling. + To run: Open python or ipython and type + "run freezingpt_slice.py". The script will prompt + you for the path to the ROMS history file, the + timestep and i-index to plot, the deepest depth + to plot, optional bounds on the colour scale, + and if 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 like. + +i_slice.py: Plot the given variable for a single i-slice (depth versus y) with + no time-averaging, spatial averaging, interpolation, or velocity + rotation. This helps to show where there are advective errors. + To run: Open python or ipython and type "run i_slice.py". The + script will prompt you for the path to the ROMS history + file, the variable name, the timestep and i-index to plot, + the deepest depth to plot, optional bounds on the colour + scale, and if 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 like. + ***ANIMATIONS*** diff --git a/freezingpt_slice.py b/freezingpt_slice.py new file mode 100644 index 0000000..e313fea --- /dev/null +++ b/freezingpt_slice.py @@ -0,0 +1,182 @@ +from netCDF4 import Dataset +from numpy import * +from matplotlib.pyplot import * +from calc_z import * + +# Plot the difference from the freezing temperature for a single i-slice +# (depth versus y) with no time-averaging, spatial averaging, or interpolation. +# This helps to show where there is spurious supercooling. +# Input: +# file_path = path to ROMS history file (NOT averages) +# tstep = timestep to plot (starting at 1) +# i_val = i-index to plot (starting at 1) +# depth_min = deepest depth to plot (negative, in metres) +# colour_bounds = optional bounds on colour scale, stored as an array of size +# 2 with the lower bound first. If colour_bounds = None, then +# determine colour scale bounds automatically. +# 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 freezingpt_slice (file_path, tstep, i_val, depth_min, colour_bounds=None, save=False, fig_name=None): + + # Grid parameters + theta_s = 0.9 + theta_b = 4.0 + hc = 40 + N = 31 + + # Read temperature, salinity, and grid variables + id = Dataset(file_path, 'r') + temp = id.variables['temp'][tstep-1,:,:-15,i_val-1] + salt = id.variables['salt'][tstep-1,:,:-15,i_val-1] + h = id.variables['h'][:-15,:] + zice = id.variables['zice'][:-15,:] + # Sea surface height is time-dependent + zeta = id.variables['zeta'][tstep-1,:-15,:] + lon_2d = id.variables['lon_rho'][:-15,:] + lat_2d = id.variables['lat_rho'][:-15,:] + id.close() + + # Calculate freezing point as seen by supercooling code + tfr = -0.054*salt + # Calculate difference from freezing point + deltat = temp - tfr + + # 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) + # Select depth and latitude at the given i-index + z = z_3d[:,:,i_val-1] + lat = tile(lat_2d[:,i_val-1], (N,1)) + + # Determine colour bounds + if colour_bounds is not None: + # Specified by user + scale_min = colour_bounds[0] + scale_max = colour_bounds[1] + if scale_min == -scale_max: + # Centered on zero; use a red-yellow-blue colour scale + colour_map = 'RdYlBu_r' + else: + # Use a rainbow colour scale + colour_map = 'jet' + else: + # Determine automatically + scale_min = amin(deltat) + scale_max = amax(deltat) + colour_map = 'jet' + + # Plot (pcolor not contour to show what each individual cell is doing) + fig = figure(figsize=(18,6)) + pcolor(lat, z, deltat, vmin=scale_min, vmax=scale_max, cmap=colour_map) + colorbar() + title('Difference from freezing point (K) at i=' + str(i_val)) + xlabel('Latitude') + ylabel('Depth (m)') + + # Determine bounds on latitude + data_sum = sum(deltat, axis=0) + edges = ma.flatnotmasked_edges(data_sum) + j_min = edges[0] + j_max = edges[1] + if j_min == 0: + # There are ocean points right to the southern boundary + # Don't do anything special + lat_min = min(lat[:,j_min]) + else: + # There is land everywhere at the southern boundary + # Show the last 2 degrees of this land mask + lat_min = min(lat[:,j_min]) - 2 + if j_max == size(data_sum) - 1: + # There are ocean points right to the northern boundary + # Don't do anything special + lat_max = max(lat[:,j_max]) + else: + # There is land everywhere at the northern boundary + # Show the first 2 degrees of this land mask + lat_max = max(lat[:,j_max]) + 2 + lat_max = -65 + xlim([lat_min, lat_max]) + ylim([depth_min, 0]) + + # Finished + if save: + fig.savefig(fig_name) + else: + fig.show() + + +# Command-line interface +if __name__ == "__main__": + + file_path = raw_input("Path to ocean history file: ") + tstep = int(raw_input("Timestep number (starting at 1): ")) + i_val = int(raw_input("i-index to plot (1 to 1443): ")) + depth_min = -1*float(raw_input("Deepest depth to plot (positive, metres): ")) + # Get colour bounds if necessary + 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 in window (d)? ") + if action == 's': + save = True + fig_name = raw_input("File name for figure: ") + elif action == 'd': + save = False + fig_name = None + + # Make the plot + freezingpt_slice(file_path, tstep, i_val, depth_min, colour_bounds, save, fig_name) + + # Repeat until the user wants to exit + while True: + repeat = raw_input("Make another plot (y/n)? ") + if repeat == 'y': + while True: + # Ask for changes to the input parameters; repeat until the user + # is finished + changes = raw_input("Enter a parameter to change: (1) file path, (2) timestep number, (3) i-index, (4) deepest depth, (5) colour bounds, (6) save/display; or enter to continue: ") + if len(changes) == 0: + # No more changes to parameters + break + else: + if int(changes) == 1: + # New file path + file_path = raw_input("Path to ocean history file: ") + elif int(changes) == 2: + # New timestep + tstep = int(raw_input("Timestep number (starting at 1): ")) + elif int(changes) == 3: + # New i-index + i_val = int(raw_input("i-index to plot (1 to 1443): ")) + elif int(changes) == 4: + # New depth bound + depth_min = -1*float(raw_input("Deepest depth to plot (positive, metres): ")) + elif int(changes) == 5: + # Get colour bounds if necessary + 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) == 6: + # Change from display to save, or vice versa + save = not save + if save: + # Get file name for figure + fig_name = raw_input("File name for figure: ") + # Make the plot + freezingpt_slice(file_path, tstep, i_val, depth_min, colour_bounds, save, fig_name) + else: + break + + + + + + + diff --git a/i_slice.py b/i_slice.py new file mode 100644 index 0000000..a04ae65 --- /dev/null +++ b/i_slice.py @@ -0,0 +1,181 @@ +from netCDF4 import Dataset +from numpy import * +from matplotlib.pyplot import * +from calc_z import * + +# Plot the given variable for a single i-slice (depth versus y) with no +# time-averaging, spatial averaging, interpolation, or velocity rotation. +# This helps to show where there are advective errors. +# Input: +# file_path = path to ROMS history file (NOT averages) +# var_name = variable name to plot +# tstep = timestep to plot (starting at 1) +# i_val = i-index to plot (starting at 1) +# depth_min = deepest depth to plot (negative, in metres) +# colour_bounds = optional bounds on colour scale, stored as an array of size +# 2 with the lower bound first. If colour_bounds = None, then +# determine colour scale bounds automatically. +# 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 i_slice (file_path, var_name, tstep, i_val, depth_min, colour_bounds=None, save=False, fig_name=None): + + # Grid parameters + theta_s = 0.9 + theta_b = 4.0 + hc = 40 + N = 31 + + # Read data and grid variables + id = Dataset(file_path, 'r') + data = id.variables[var_name][tstep-1,:,:-15,i_val-1] + h = id.variables['h'][:-15,:] + zice = id.variables['zice'][:-15,:] + # Sea surface height is time-dependent + zeta = id.variables['zeta'][tstep-1,:-15,:] + lon_2d = id.variables['lon_rho'][:-15,:] + lat_2d = id.variables['lat_rho'][:-15,:] + id.close() + + # Get a 3D array of z-coordinates; sc_r and Cs_r are unused in this script + z_3d, sc_r, Cs_r = calc_z(h, zice, theta_s, theta_b, hc, N, zeta) + # Select depth and latitude at the given i-index + z = z_3d[:,:,i_val-1] + lat = tile(lat_2d[:,i_val-1], (N,1)) + + # Determine colour bounds + if colour_bounds is not None: + # Specified by user + scale_min = colour_bounds[0] + scale_max = colour_bounds[1] + if scale_min == -scale_max: + # Centered on zero; use a red-yellow-blue colour scale + colour_map = 'RdYlBu_r' + else: + # Use a rainbow colour scale + colour_map = 'jet' + else: + # Determine automatically + scale_min = amin(data) + scale_max = amax(data) + colour_map = 'jet' + + # Plot (pcolor not contour to show what each individual cell is doing) + fig = figure(figsize=(18,6)) + pcolor(lat, z, data, vmin=scale_min, vmax=scale_max, cmap=colour_map) + colorbar() + title(var_name + ' at i=' + str(i_val)) + xlabel('Latitude') + ylabel('Depth (m)') + + # Determine bounds on latitude + data_sum = sum(data, axis=0) + edges = ma.flatnotmasked_edges(data_sum) + j_min = edges[0] + j_max = edges[1] + if j_min == 0: + # There are ocean points right to the southern boundary + # Don't do anything special + lat_min = min(lat[:,j_min]) + else: + # There is land everywhere at the southern boundary + # Show the last 2 degrees of this land mask + lat_min = min(lat[:,j_min]) - 2 + if j_max == size(data_sum) - 1: + # There are ocean points right to the northern boundary + # Don't do anything special + lat_max = max(lat[:,j_max]) + else: + # There is land everywhere at the northern boundary + # Show the first 2 degrees of this land mask + lat_max = max(lat[:,j_max]) + 2 + lat_max = -65 + xlim([lat_min, lat_max]) + ylim([depth_min, 0]) + + # Finished + if save: + fig.savefig(fig_name) + else: + fig.show() + + +# Command-line interface +if __name__ == "__main__": + + file_path = raw_input("Path to ocean history file: ") + var_name = raw_input("Variable name: ") + tstep = int(raw_input("Timestep number (starting at 1): ")) + i_val = int(raw_input("i-index to plot (1 to 1443): ")) + depth_min = -1*float(raw_input("Deepest depth to plot (positive, metres): ")) + # Get colour bounds if necessary + 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 in window (d)? ") + if action == 's': + save = True + fig_name = raw_input("File name for figure: ") + elif action == 'd': + save = False + fig_name = None + + # Make the plot + i_slice(file_path, var_name, tstep, i_val, depth_min, colour_bounds, save, fig_name) + + # Repeat until the user wants to exit + while True: + repeat = raw_input("Make another plot (y/n)? ") + if repeat == 'y': + while True: + # Ask for changes to the input parameters; repeat until the user + # is finished + changes = raw_input("Enter a parameter to change: (1) file path, (2) variable name, (3) timestep number, (4) i-index, (5) deepest depth, (6) colour bounds, (7) save/display; or enter to continue: ") + if len(changes) == 0: + # No more changes to parameters + break + else: + if int(changes) == 1: + # New file path + file_path = raw_input("Path to ocean history file: ") + elif int(changes) == 2: + # New variable name + var_name = raw_input("Variable name: ") + elif int(changes) == 3: + # New timestep + tstep = int(raw_input("Timestep number (starting at 1): ")) + elif int(changes) == 4: + # New i-index + i_val = int(raw_input("i-index to plot (1 to 1443): ")) + elif int(changes) == 5: + # New depth bound + depth_min = -1*float(raw_input("Deepest depth to plot (positive, metres): ")) + elif int(changes) == 6: + # Get colour bounds if necessary + 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) == 7: + # Change from display to save, or vice versa + save = not save + if save: + # Get file name for figure + fig_name = raw_input("File name for figure: ") + # Make the plot + i_slice(file_path, var_name, tstep, i_val, depth_min, colour_bounds, save, fig_name) + else: + break + + + + + + + diff --git a/ismr_map.py b/ismr_map.py index 1c4de7e..c24677c 100644 --- a/ismr_map.py +++ b/ismr_map.py @@ -21,8 +21,8 @@ def ismr_map (grid_path, log_path, save=False, fig_name=None): lon_max = [-59.33, -60, -66.67, -28.33, -88.83, -99.17, -103.33, -111.5, -114.33, -140, -145, 146.62, 123.33, 102.5, 89.17, 75, 37.67, 33.33, 16.17, 12.88, 7.6, -10.33, -146.67, 180] lat_min = [-73.03, -69.35, -74.17, -83.5, -73.28, -75.5, -75.5, -75.33, -74.9, -76.42, -78, -67.83, -67.17, -66.67, -67.83, -73.67, -69.83, -71.67, -70.5, -70.75, -71.83, -76.33, -85, -84.5] lat_max = [-69.37, -66.13, -69.5, -74.67, -71.67, -74.17, -74.67, -73.67, -73, -75.17, -76.41, -66.67, -66.5, -64.83, -66.17, -68.33, -68.67, -68.33, -69.33, -69.83, -69.33, -71.5, -77.77, -77] - # Area of each ice shelf in m^2 (printed to screen during massloss.py, - # update if the grid changes) + # Area of each ice shelf in m^2 (printed to screen during + # timeseries_massloss.py, update if the grid changes) area = [12754001970.4, 52008964915.9, 47287926558.6, 353435138251.0, 31290573250.5, 5162051654.52, 3990382861.08, 4680996769.75, 32446806852.2, 7694313052.38, 13537287121.0, 4918446447.87, 6482036686.01, 30521756982.6, 15158334399.6, 64359735004.9, 4575785549.65, 45327465354.5, 8110511960.62, 7088165282.99, 54898163328.1, 68006982027.4, 429252991746.0] # Observed melt rate (Rignot 2013) and uncertainty for each ice shelf, in Gt/y obs_ismr = [0.1, 0.4, 3.1, 0.3, 1.7, 16.2, 17.7, 7.8, 4.3, 0.6, 1.5, 1.4, 7.7, 2.8, 1.7, 0.6, -0.4, 0.4, 0.7, 0.5, 0.5, 0.1, 0.1] @@ -56,7 +56,7 @@ def ismr_map (grid_path, log_path, save=False, fig_name=None): # Skip the values for the entire continent for line in f: try: - pass + tmp = float(line) except(ValueError): break # Set up array for melt rate values at each ice shelf @@ -130,7 +130,7 @@ def ismr_map (grid_path, log_path, save=False, fig_name=None): error[region] = error_tmp # Edit zice so tiny ice shelves won't be contoured - zice[error.mask] = 0.0 + #zice[error.mask] = 0.0 # Convert grid to spherical coordinates x = -(lat+90)*cos(lon*deg2rad+pi/2) diff --git a/massloss_map.py b/massloss_map.py index e4f2935..19e55d3 100644 --- a/massloss_map.py +++ b/massloss_map.py @@ -17,8 +17,8 @@ def massloss_map (grid_path, log_path, save=False, fig_name=None): # These depend on the source geometry, in this case RTopo 1.05 # Note there is one extra index at the end of each array; this is because # the Ross region crosses the line 180W and therefore is split into two - lon_min = [-62.67, -65.5, -79.17, -85, -104.17, -102.5, -108.33, -114.5, -135.67, -149.17, -155, 144, 115, 94.17, 80.83, 65, 33.83, 19, 12.9, 9.33, -10.05, -28.33, -180, 158.33] - lon_max = [-59.33, -60, -66.67, -28.33, -88.83, -99.17, -103.33, -111.5, -114.33, -140, -145, 146.62, 123.33, 102.5, 89.17, 75, 37.67, 33.33, 16.17, 12.88, 7.6, -10.33, -146.67, 180] + lon_min = [-62.67, -65.5, -79.17, -85, -104.17, -102.5, -108.33, -114.5, -135.67, -149.17, -155, 144, 115, 94.17, 80.83, 65, 33.83, 19, 12.9, 9.33, -10.05, -28.33, -181, 158.33] + lon_max = [-59.33, -60, -66.67, -28.33, -88.83, -99.17, -103.33, -111.5, -114.33, -140, -145, 146.62, 123.33, 102.5, 89.17, 75, 37.67, 33.33, 16.17, 12.88, 7.6, -10.33, -146.67, 181] lat_min = [-73.03, -69.35, -74.17, -83.5, -73.28, -75.5, -75.5, -75.33, -74.9, -76.42, -78, -67.83, -67.17, -66.67, -67.83, -73.67, -69.83, -71.67, -70.5, -70.75, -71.83, -76.33, -85, -84.5] lat_max = [-69.37, -66.13, -69.5, -74.67, -71.67, -74.17, -74.67, -73.67, -73, -75.17, -76.41, -66.67, -66.5, -64.83, -66.17, -68.33, -68.67, -68.33, -69.33, -69.83, -69.33, -71.5, -77.77, -77] # Observed mass loss (Rignot 2013) and uncertainty for each ice shelf, in Gt/y @@ -51,7 +51,7 @@ def massloss_map (grid_path, log_path, save=False, fig_name=None): # Skip the values for the entire continent for line in f: try: - pass + tmp = float(line) except(ValueError): break # Set up array for mass loss values at each ice shelf @@ -65,7 +65,7 @@ def massloss_map (grid_path, log_path, save=False, fig_name=None): massloss_ts[index, t] = float(line) t += 1 except(ValueError): - # Reached the header for the next variable + # Reached the header for the next ice shelf break index +=1 @@ -124,7 +124,7 @@ def massloss_map (grid_path, log_path, save=False, fig_name=None): error[region] = error_tmp # Edit zice so tiny ice shelves won't be contoured - zice[error.mask] = 0.0 + #zice[error.mask] = 0.0 # Convert grid to spherical coordinates x = -(lat+90)*cos(lon*deg2rad+pi/2) diff --git a/mmm_atmos_netcdf.py b/mmm_atmos_netcdf.py deleted file mode 100644 index 9cd4cc9..0000000 --- a/mmm_atmos_netcdf.py +++ /dev/null @@ -1,100 +0,0 @@ -from cmip5_paths import * -from numpy import * -from netCDF4 import Dataset - -# Calculate the multi-model mean of atmospheric climatology files created using -# cmip5_atmos_climatology_netcdf.py. -def mmm_atmos_netcdf (): - - # Directory containing CMIP5 NetCDF climatology atmosphere files, - # interpolated to the ERA-Interim grid (created using - # cmip5_atmos_climatology_netcdf.py) - directory = '/short/y99/kaa561/CMIP5_forcing/atmos/' - # Path to ERA-Interim file (created using eraint_climatology_netcdf.py) - eraint_file = directory + 'ERA-Interim.nc' - # Path to multi-model mean output file - output_file = directory + 'MMM.nc' - # Variable names in NetCDF files - var_names = ['Pair', 'Tair', 'Hair', 'cloud', 'Uwind', 'Vwind', 'precip', 'snow', 'evap', 'swrad', 'lwrad'] - # Corresponding units - var_units = ['kPa', 'degC', '1', '%', 'm/s', 'm/s', '10^6 kg/m^2/s', '10^6 kg/m^2/s', '10^6 kg/m^2/s', 'W/m^2', 'W/m^2'] - - # Read ERA-Interim grid - id = Dataset(eraint_file, 'r') - lon = id.variables['longitude'][:] - lat = id.variables['latitude'][:] - id.close() - - # Set up output file - print 'Setting up ' + output_file - out_id = Dataset(output_file, 'w') - # Define dimensions - out_id.createDimension('longitude', size(lon)) - out_id.createDimension('latitude', size(lat)) - out_id.createDimension('time', 12) - # Define dimension variables and fill with axes - out_id.createVariable('longitude', 'f8', ('longitude')) - out_id.variables['longitude'].units = 'degrees' - out_id.variables['longitude'][:] = lon - out_id.createVariable('latitude', 'f8', ('latitude')) - out_id.variables['latitude'].units = 'degrees' - out_id.variables['latitude'][:] = lat - out_id.createVariable('time', 'f8', ('time')) - out_id.variables['time'].units = 'month' - out_id.variables['time'][:] = arange(1, 12+1) - - # Get a list of CMIP5 Model objects - models = build_model_list() - # Build a corresponding list of model names - model_names = [] - for model in models: - model_names.append(model.name) - - # Loop over variables - for i in range(len(var_names)): - var = var_names[i] - print 'Variable ' + var - - num_models = 0 # Number of models in the multi-model mean - multi_model_mean = None - # Loop over models - for model_name in model_names: - # Read model data - id = Dataset(directory + model_name + '.nc', 'r') - model_data = id.variables[var][:,:,:] - id.close() - - # Check for missing data - try: - mask = model_data.mask - except(AttributeError): - # There is no mask; set it to False - mask = False - if all(mask): - # Everything is masked; data is missing for this variable - pass - else: - # Add to multi-model mean and increment num_models - if multi_model_mean is None: - multi_model_mean = model_data[:,:,:] - else: - multi_model_mean[:,:,:] += model_data[:,:,:] - num_models += 1 - - # Divide multi_model_mean (currently just a sum) by num_models - multi_model_mean /= num_models - # Define variable and fill with data - out_id.createVariable(var_names[i], 'f8', ('time', 'latitude', 'longitude')) - out_id.variables[var_names[i]].units = var_units[i] - out_id.variables[var_names[i]][:,:,:] = multi_model_mean - out_id.close() - - -# Command-line interface -if __name__ == "__main__": - - mmm_atmos_netcdf() - - - - diff --git a/mmm_ocean_netcdf.py b/mmm_ocean_netcdf.py deleted file mode 100644 index 7424f45..0000000 --- a/mmm_ocean_netcdf.py +++ /dev/null @@ -1,94 +0,0 @@ -from cmip5_paths import * -from numpy import * -from netCDF4 import Dataset - -# Calculate the multi-model mean of ocean climatology files created using -# cmip5_ocean_climatology_netcdf.py. -def mmm_ocean_netcdf (): - - # Directory containing CMIP5 NetCDF climatology ocean files, interpolated - # to the ECCO2 grid at the northern boundary of ROMS (created using - # cmip5_ocean_climatology_netcdf.py) - directory = '/short/y99/kaa561/CMIP5_forcing/ocean/' - # Path to ECCO2 file (created using cmip5_ocean_climatology_netcdf.py) - ecco2_file = directory + 'ECCO2.nc' - # Path to multi-model mean output file - output_file = directory + 'MMM.nc' - # Variable names in NetCDF files - var_names = ['temp', 'salt', 'u', 'v'] - # Corresponding units - var_units = ['degC', 'psu', 'm/s', 'm/s'] - - # Read ECCO2 grid - id = Dataset(ecco2_file, 'r') - lon = id.variables['longitude'][:] - depth = id.variables['depth'][:] - id.close() - - # Set up output file - print 'Setting up ' + output_file - out_id = Dataset(output_file, 'w') - # Define dimensions - out_id.createDimension('longitude', size(lon)) - out_id.createDimension('depth', size(depth)) - out_id.createDimension('time', 12) - # Define dimension variables and fill with axes - out_id.createVariable('longitude', 'f8', ('longitude')) - out_id.variables['longitude'].units = 'degrees' - out_id.variables['longitude'][:] = lon - out_id.createVariable('depth', 'f8', ('depth')) - out_id.variables['depth'].units = 'metres' - out_id.variables['depth'][:] = depth - out_id.createVariable('time', 'f8', ('time')) - out_id.variables['time'].units = 'month' - out_id.variables['time'][:] = arange(1, 12+1) - - # Get a list of CMIP5 Model objects - models = build_model_list() - # Build a corresponding list of model names - model_names = [] - for model in models: - model_names.append(model.name) - - # Loop over variables - for i in range(len(var_names)): - var = var_names[i] - print 'Variable ' + var - - num_models = 0 # Number of models in the multi-model mean - multi_model_mean = None - # Loop over models - for model_name in model_names: - # Read model data - id = Dataset(directory + model_name + '.nc', 'r') - model_data = id.variables[var][:,:,:] - id.close() - - if all(model_data.mask): - # Everything is masked; data is missing for this variable - pass - else: - # Add to multi-model mean and increment num_models - if multi_model_mean is None: - multi_model_mean = model_data[:,:,:] - else: - multi_model_mean[:,:,:] += model_data[:,:,:] - num_models += 1 - - # Divide multi_model_mean (currently just a sum) by num_models - multi_model_mean /= num_models - # Define variable and fill with data - out_id.createVariable(var_names[i], 'f8', ('time', 'depth', 'longitude')) - out_id.variables[var_names[i]].units = var_units[i] - out_id.variables[var_names[i]][:,:,:] = multi_model_mean - out_id.close() - - -# Command-line interface -if __name__ == "__main__": - - mmm_ocean_netcdf() - - - - diff --git a/rignot_data b/rignot_data deleted file mode 100644 index 7a064b1..0000000 --- a/rignot_data +++ /dev/null @@ -1,25 +0,0 @@ -Name Area (km2) Mass loss (Gt/y) Melt rate (m/y) lon-bounds lat-bounds ROMS area (km^2) Area bias (%) -Larsen D 22,548 1.4 +/- 14 0.1 +/- 0.6 -62.67:-59.33 -73.03:-69.37 12,754 -43 -Larsen C 46,465 20.7 +/- 67 0.4 +/- 1 -65.5:-60 -69.35:-66.13 52,009 12 -Wilkins & George VI & Stange 44,327 135.4 +/- 40 3.1 +/- 0.8 -79.17:-66.67 -74.17:-69.5 47,287 7 -Ronne-Filchner 443,140 155.4 +/- 45 0.3 +/- 0.1 -85:-28.33 -83.5:-74.67 353,435 -20 -Abbot 29,688 51.8 +/- 19 1.7 +/- 0.6 -104.17:-88.83 -73.28:-71.67 31,291 5 -Pine Island 6,249 101.2 +/- 8 16.2 +/- 1 -102.5:-99.17 -75.5:-74.17 5,162 -17 -Thwaites 5,499 97.5 +/- 7 17.7 +/- 1 -108.33:-103.33 -75.5:-74.67 3,990 -27 -Dotson 5,803 45.2 +/- 4 7.8 +/- 0.6 -114.5:-111.5 -75.33:-73.67 4,681 -19 -Getz 34,018 144.9 +/- 14 4.3 +/- 0.4 -135.67:-114.33 -74.9:-73 32,447 -5 -Nickerson 6,495 4.2 +/- 2 0.6 +/- 0.3 -149.17:-140 -76.42:-75.17 7,694 18 -Sulzberger 12,333 18.2 +/- 3 1.5 +/- 0.3 -155:-145 -78:-76.41 13,537 10 -Ross 500,809 47.7 +/- 34 0.1 +/- 0.1 -180:-146.67,158.33:180 -85:-77.77,-84.5:-77 429,253 -14 -Mertz 5,522 7.9 +/- 3 1.4 +/- 0.6 144:146.62 -67.83:-66.67 4,918 -11 -Totten & Moscow Uni 11,830 90.6 +/- 8 7.7 +/- 0.7 115:123.33 -67.17:-66.5 6,482 -45 -Shackleton 26,080 72.6 +/- 15 2.8 +/- 0.6 94.17:102.5 -66.67:-64.83 30,522 17 -West 15,666 27.2 +/- 10 1.7 +/- 0.7 80.83:89.17 -67.83:-66.17 15,158 -3 -Amery 60,654 35.5 +/- 23 0.6 +/- 0.4 65:75 -73.67:-68.33 64,360 6 -Prince Harald 5,392 -2 +/- 3 -0.4 +/- 0.6 33.83:37.67 -69.83:-68.67 4,576 -15 -Baudouin & Borchgrevink 54,532 21.6 +/- 18 0.4 +/- 0.4 19:33.33 -71.67:-68.33 45,327 -17 -Lazarev 8,519 6.3 +/- 2 0.7 +/- 0.2 12.9:16.17 -70.5:-69.33 8,111 -5 -Nivl 7,285 3.9 +/- 2 0.5 +/- 0.2 9.33:12.88 -70.75:-69.83 7,088 -3 -Fimbul & Jelbart & Ekstrom 58,559 26.8 +/- 14 0.5 +/- 0.2 -10.05:7.6 -71.83:-69.33 55,763 -5 -Brunt & Riiser-Larsen 80,344 9.7 +/- 16 0.1 +/- 0.2 -28.33:-10.33 -76.33:-71.5 68,007 -15 -