Added RMS tidal velocity code; added FESOM high-res to mip_grid_res.py

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

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);

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,15 +109,16 @@ 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
         fig_name = raw_input("File name for figure: ")
     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)