From 01a12fdc06335c3908252cd0c25a19ef5b9bf416 Mon Sep 17 00:00:00 2001 From: Kaitlin Naughten Date: Fri, 7 Jul 2017 16:50:44 +1000 Subject: [PATCH] Added RMS tidal velocity code; added FESOM high-res to mip_grid_res.py --- file_guide.txt | 18 ++++++++---- get_rms_tide_vel.m | 57 ++++++++++++++++++++++++++++++++++++ mip_grid_res.py | 73 +++++++++++++++++++++++++++++----------------- 3 files changed, 117 insertions(+), 31 deletions(-) create mode 100644 get_rms_tide_vel.m diff --git a/file_guide.txt b/file_guide.txt index 797b861..d2dac55 100644 --- a/file_guide.txt +++ b/file_guide.txt @@ -300,6 +300,14 @@ romscice_flux_correction.py: Build a monthly climatology of the extra surface the concatenated ROMS averages file and the desired output forcing file. +get_rms_tide_vel.m: Matlab script (the horror!) which calls the TMD/CATS2008a + tidal model: + https://www.esr.org/polar_tide_models/Model_CATS2008a.html + to get the RMS tidal velocity (sum of 10 constituents) for + each point on the ROMS grid. Outputs to a NetCDF file, which + can be used as a forcing file in MetROMS if you activate + ICESHELF_RMS_TIDES. + ***POST-PROCESSING DIAGNOSTIC FILES*** @@ -1124,18 +1132,18 @@ gadv_frazil_bathy.py: Plot a section of the coastline from 25-45E. On the left, ***FIGURES FOR FESOM INTERCOMPARISON*** -mip_grid_res.py: Make a 2x1 plot showing horizontal grid resolution (square +mip_grid_res.py: Make a 3x1 plot showing horizontal grid resolution (square root of the area of each rectanuglar grid cell or triangular - element) in MetROMS (left) and FESOM (right), zoomed into the - Antarctic continental shelf. + element) in MetROMS, FESOM low-res, and FESOM high-res, zoomed + into the Antarctic continental shelf. To run: First clone my "fesomtools" GitHub repository and replace '/short/y99/kaa561/fesomtools' (near the top of the file) with the path to the cloned repository on your system. Then open python or ipython and type "run mip_grid_res.py". The script will prompt you for the path to the ROMS grid file and FESOM mesh - directory, and whether you want to save the figure (and - if so, what filename) or display it on the screen. + directories, and whether you want to save the figure + (and if so, what filename) or display it on the screen. mip_timeseries.py: Plot MetROMS and FESOM timeseries together, for Drake Passage transport, total Antarctic sea ice area, volume, and diff --git a/get_rms_tide_vel.m b/get_rms_tide_vel.m new file mode 100644 index 0000000..b807661 --- /dev/null +++ b/get_rms_tide_vel.m @@ -0,0 +1,57 @@ +% Read ROMS grid +id = netcdf.open('circ30S_quarterdegree.nc','NOWRITE'); +lat_id = netcdf.inqVarID(id, 'lat_rho'); +lat = netcdf.getVar(id, lat_id); +lon_id = netcdf.inqVarID(id, 'lon_rho'); +lon = netcdf.getVar(id, lon_id); +h_id = netcdf.inqVarID(id, 'h'); +h = netcdf.getVar(id, h_id); +zice_id = netcdf.inqVarID(id, 'zice'); +zice = netcdf.getVar(id, zice_id); +netcdf.close(id); + +% Make sure longitude is in the range (-180, 180) +index = lon > 180; +lon(index) = lon(index)-360; +num_lon = size(lon,1); +num_lat = size(lon,2); + +% Get water column thickness +wct = h + zice; + +% Call TMD for the amplitudes of tidal transports U and V in m^2/s +% for each tidal constituent +[amp_U,~,~,~] = tmd_extract_HC('DATA/Model_CATS2008a_opt', lat, lon, 'U'); +[amp_V,~,~,~] = tmd_extract_HC('DATA/Model_CATS2008a_opt', lat, lon, 'V'); +% Get directionless amplitude with sqrt(amp_U^2 + amp_V^2) +% Divide by sqrt(2) for RMS of sinusoid +rms_components = sqrt(amp_U.^2 + amp_V.^2)/sqrt(2); +% Sum the tidal consitutents and divide by water column thickness for +% result in m/s +rms_tide_vel_tmp = squeeze(sum(rms_components, 1))./wct; +% Add a time axis of size 1 +rms_tide_vel = zeros([1, size(rms_tide_vel_tmp,1), size(rms_tide_vel_tmp,2)]); +rms_tide_vel(1,:,:) = rms_tide_vel_tmp; + +% Remove NaNs +index = isnan(rms_tide_vel); +rms_tide_vel(index) = 0.0; +% Set to zero outside ice shelf cavities +index = zice==0; +rms_tide_vel(index) = 0; + +% Write to file +id = netcdf.create('rms_tides.nc', 'CLOBBER'); +time_dim = netcdf.defDim(id, 'rms_time', netcdf.getConstant('NC_UNLIMITED')); +eta_dim = netcdf.defDim(id, 'eta_rho', num_lat); +xi_dim = netcdf.defDim(id, 'xi_rho', num_lon); +time_id = netcdf.defVar(id, 'rms_time', 'double', time_dim); +netcdf.putAtt(id, time_id, 'calendar', 'none'); +tide_id = netcdf.defVar(id, 'tide_RMSvel', 'NC_DOUBLE', [xi_dim, eta_dim, time_dim]); +netcdf.putAtt(id, tide_id, 'units', 'm/s'); +netcdf.endDef(id); +netcdf.putVar(id, time_id, 0, 1, [0.]); +start = [0, 0, 0]; +count = [num_lon, num_lat, 1]; +netcdf.putVar(id, tide_id, start, count, rms_tide_vel); +netcdf.close(id); diff --git a/mip_grid_res.py b/mip_grid_res.py index 056f360..7adccf5 100644 --- a/mip_grid_res.py +++ b/mip_grid_res.py @@ -8,16 +8,17 @@ sys.path.insert(0, '/short/y99/kaa561/fesomtools') from patches import * -# Make a 2x1 plot showing horizontal grid resolution (square root of the area -# of each rectanuglar grid cell or triangular element) in MetROMS (left) and -# FESOM (right), zoomed into the Antarctic continental shelf. +# Make a 3x1 plot showing horizontal grid resolution (square root of the area +# of each rectanuglar grid cell or triangular element) in MetROMS, FESOM +# low-res, and FESOM high-res, zoomed into the Antarctic continental shelf. # Input: # roms_grid_file = path to ROMS grid file -# fesom_mesh_dir = path to FESOM mesh directory +# fesom_mesh_low, fesom_mesh_high = paths to FESOM low-res and high-res mesh +# directories # save = optional boolean indicating to save the figure, rather than display it # on screen # fig_name = if save=True, filename for figure -def mip_grid_res (roms_grid_file, fesom_mesh_dir, save=False, fig_name=None): +def mip_grid_res (roms_grid_file, fesom_mesh_low, fesom_mesh_high, save=False, fig_name=None): # Spatial bounds on plot lat_max = -63 + 90 @@ -45,38 +46,57 @@ def mip_grid_res (roms_grid_file, fesom_mesh_dir, save=False, fig_name=None): roms_x = -(roms_lat+90)*cos(roms_lon*deg2rad+pi/2) roms_y = (roms_lat+90)*sin(roms_lon*deg2rad+pi/2) - print 'Processing FESOM' + print 'Processing FESOM low-res' # Build triangular patches for each element - elements, patches = make_patches(fesom_mesh_dir, circumpolar) + elements_low, patches_low = make_patches(fesom_mesh_low, circumpolar) # Calculate the resolution at each element - fesom_res = [] - for elm in elements: - fesom_res.append(sqrt(elm.area())*1e-3) + fesom_res_low = [] + for elm in elements_low: + fesom_res_low.append(sqrt(elm.area())*1e-3) + + print 'Processing FESOM high-res' + # Build triangular patches for each element + elements_high, patches_high = make_patches(fesom_mesh_high, circumpolar) + # Calculate the resolution at each element + fesom_res_high = [] + for elm in elements_high: + fesom_res_high.append(sqrt(elm.area())*1e-3) print 'Plotting' - fig = figure(figsize=(20,9)) + fig = figure(figsize=(27,9)) # ROMS - ax1 = fig.add_subplot(1,2,1, aspect='equal') + ax1 = fig.add_subplot(1,3,1, aspect='equal') pcolor(roms_x, roms_y, roms_res, vmin=limits[0], vmax=limits[1], cmap='jet') xlim([-lat_max, lat_max]) ylim([-lat_max, lat_max]) axis('off') - title('MetROMS', fontsize=24) - # FESOM - ax2 = fig.add_subplot(1,2,2, aspect='equal') - img = PatchCollection(patches, cmap='jet') - img.set_array(array(fesom_res)) - img.set_clim(vmin=limits[0], vmax=limits[1]) - img.set_edgecolor('face') - ax2.add_collection(img) + title('MetROMS', fontsize=28) + # FESOM low-res + ax2 = fig.add_subplot(1,3,2, aspect='equal') + img_low = PatchCollection(patches_low, cmap='jet') + img_low.set_array(array(fesom_res_low)) + img_low.set_clim(vmin=limits[0], vmax=limits[1]) + img_low.set_edgecolor('face') + ax2.add_collection(img_low) + xlim([-lat_max, lat_max]) + ylim([-lat_max, lat_max]) + axis('off') + title('FESOM low-res', fontsize=28) + # FESOM high-res + ax3 = fig.add_subplot(1,3,3, aspect='equal') + img_high = PatchCollection(patches_high, cmap='jet') + img_high.set_array(array(fesom_res_high)) + img_high.set_clim(vmin=limits[0], vmax=limits[1]) + img_high.set_edgecolor('face') + ax3.add_collection(img_high) xlim([-lat_max, lat_max]) ylim([-lat_max, lat_max]) axis('off') - title('FESOM', fontsize=24) + title('FESOM high-res', fontsize=28) cbaxes = fig.add_axes([0.92, 0.2, 0.01, 0.6]) - cbar = colorbar(img, cax=cbaxes, extend='max', ticks=arange(limits[0], limits[1]+5, 5)) - cbar.ax.tick_params(labelsize=20) - suptitle('Horizontal grid resolution (km)', fontsize=30) + cbar = colorbar(img_high, cax=cbaxes, extend='max', ticks=arange(limits[0], limits[1]+5, 5)) + cbar.ax.tick_params(labelsize=24) + suptitle('Horizontal grid resolution (km)', fontsize=36) subplots_adjust(wspace=0.05) if save: @@ -89,7 +109,8 @@ def mip_grid_res (roms_grid_file, fesom_mesh_dir, save=False, fig_name=None): if __name__ == "__main__": roms_grid_file = raw_input("Path to ROMS grid file: ") - fesom_mesh_dir = raw_input("Path to FESOM mesh directory: ") + fesom_mesh_low = raw_input("Path to FESOM low-res mesh directory: ") + fesom_mesh_high = raw_input("Path to FESOM high-res mesh directory: ") action = raw_input("Save figure (s) or display in window (d)? ") if action == 's': save = True @@ -97,7 +118,7 @@ def mip_grid_res (roms_grid_file, fesom_mesh_dir, save=False, fig_name=None): elif action == 'd': save = False fig_name = None - mip_grid_res (roms_grid_file, fesom_mesh_dir, save, fig_name) + mip_grid_res (roms_grid_file, fesom_mesh_low, fesom_mesh_high, save, fig_name)