From 7d43ead3e7e914f8f03be55024a86e5796fca9ea Mon Sep 17 00:00:00 2001 From: Kaitlin Naughten Date: Mon, 29 May 2017 16:37:37 +1000 Subject: [PATCH] Comments and documentation --- adv_mld_2.py | 6 +- adv_thickness_2.py | 5 +- calc_rx1.py | 18 +++++- file_guide.txt | 125 +++++++++++++++++++++++++++++++++------- find_isolated_points.py | 16 ++++- fix_isolated_points.py | 13 +++++ plot_layers.py | 22 +++++++ romscice_sss_nudging.py | 42 +++++++++++++- 8 files changed, 218 insertions(+), 29 deletions(-) diff --git a/adv_mld_2.py b/adv_mld_2.py index 53524cc..ac58991 100644 --- a/adv_mld_2.py +++ b/adv_mld_2.py @@ -3,9 +3,9 @@ from matplotlib.pyplot import * #import colormaps as cmaps -# Create a 3x2 plot of mixed layer depth (calculated by KPP) on 23 August (the -# sea ice area max) for each advection experiment: the absolute values for -# U3_LIM, and the anomalies from U3_LIM for the other 5 experiments. +# Create a 2x1 plot of mixed layer depth (calculated by KPP) on 23 August (the +# sea ice area max) for the U3_LIM experiment, and anomalies for the C4_LD +# experiment. def adv_mld_2 (): # Paths to simulation directories diff --git a/adv_thickness_2.py b/adv_thickness_2.py index ad1316b..fda647f 100644 --- a/adv_thickness_2.py +++ b/adv_thickness_2.py @@ -3,9 +3,8 @@ from matplotlib.pyplot import * #import colormaps as cmaps -# Create a 3x2 plot of sea ice effective thickness on 23 August (the sea ice -# area max) for each advection experiment: the absolute values for U3_LIM, and -# the anomalies from U3_LIM for the other 5 experiments. +# Create a 2x1 plot of sea ice effective thickness on 23 August (the sea ice +# area max) for the U3_LIM experiment, and anomalies for the C4_LD experiment. def adv_thickness_2 (): # Paths to simulation directories diff --git a/calc_rx1.py b/calc_rx1.py index 111cc95..b3518d7 100644 --- a/calc_rx1.py +++ b/calc_rx1.py @@ -2,8 +2,15 @@ from numpy import * from calc_z import * +# Calculate the rx1 field (measure of grid steepness) for the given grid and +# vertical stretching parameters. Save to a NetCDF file. +# Input: +# grid_path = path to ROMS grid file +# theta_s, theta_b, hc, N, Vstretching = grid stretching parameters (check *.in) +# out_file = path to desired output file def calc_rx1 (grid_path, theta_s, theta_b, hc, N, Vstretching, out_file): + # Read grid id = Dataset(grid_path, 'r') lon_2d = id.variables['lon_rho'][:,:] lat_2d = id.variables['lat_rho'][:,:] @@ -12,23 +19,31 @@ def calc_rx1 (grid_path, theta_s, theta_b, hc, N, Vstretching, out_file): mask_rho = id.variables['mask_rho'][:,:] id.close() + # Calculate 3D field of depth values z, sc_r, Cs_r = calc_z(h, zice, theta_s, theta_b, hc, N, None, Vstretching) + # Mask out land for k in range(N): tmp = z[k,:,:] tmp[mask_rho==0] = NaN z[k,:,:] = tmp + # Calculate rx1 in each dimension rx1_3d_i = abs((z[1:,1:,1:]-z[1:,1:,:-1]+z[:-1,1:,1:]-z[:-1,1:,:-1])/(z[1:,1:,1:]+z[1:,1:,:-1]-z[:-1,1:,1:]-z[:-1,1:,:-1])) rx1_3d_j = abs((z[1:,1:,1:]-z[1:,:-1,1:]+z[:-1,1:,1:]-z[:-1,:-1,1:])/(z[1:,1:,1:]+z[1:,:-1,1:]-z[:-1,1:,1:]-z[:-1,:-1,1:])) + # Save the larger of the two rx1_3d = maximum(rx1_3d_i, rx1_3d_j) + # Take maximum along depth to get a 2D field rx1_tmp = amax(rx1_3d, axis=0) + # This only worked for interior points; copy the boundary in each dimesion rx1 = zeros(shape(lon_2d)) rx1[1:,1:] = rx1_tmp rx1[0,:] = rx1[1,:] rx1[:,0] = rx1[:,1] + # Mask out NaNs rx1 = ma.masked_where(isnan(rx1), rx1) + # Write output file id = Dataset(out_file, 'w') id.createDimension('xi_rho', size(lon_2d,1)) id.createDimension('eta_rho', size(lon_2d,0)) @@ -36,7 +51,8 @@ def calc_rx1 (grid_path, theta_s, theta_b, hc, N, Vstretching, out_file): id.variables['rx1'][:,:] = rx1 id.close() - + +# Command-line interface if __name__ == "__main__": grid_path = raw_input("Path to grid file: ") diff --git a/file_guide.txt b/file_guide.txt index c24b1d3..fcd4643 100644 --- a/file_guide.txt +++ b/file_guide.txt @@ -262,6 +262,52 @@ romscice_gpcp.py: Convert one year of GPCP monthly averages for total and type "run romscice_gpcp.py". The years 1992-2005 will be processed in a loop. +find_isolated_points.py: Find all of the CICE grid points which are land (or ice + shelf) on 3 sides. Sea ice can grow in these isolated + points but cannot escape due to CICE's coastal boundary + conditions, so it gets crazy thick (like 2 km thick). + Print the indices of these points to the screen. This + script assumes a periodic boundary in the longitude + direction + To run: Edit the variables start_j, end_j for the range + of j-indices you care about (i.e. could + possibly grow sea ice). Then open python or + ipython and type "run find_isolated_points.py". + The script will prompt you for the path to your + CICE land mask (kmt) file. Run cice_grid.py + if this file doesn't exist yet. + +fix_isolated_points.py: Based on the results of find_isolated_points.py, fix + all the points which are isolated on 3 sides in the CICE + grid, by editing the ROMS grid. Set some of them to ice + shelf points (with zice the average of the correct + neighbour values), some to land points, remove the ice + shelves on others, and make other land points ocean + points (with h the average of the correct neighbour + values). After running this script, rerun cice_grid.py + to update the CICE grid. This script is grid-specific + and manually written. If you are editing it for a new + grid, make sure you take into account the ghost cells + in ROMS and CICE: point (i,j) in CICE is point (i+1,j+1) + in ROMS. + To run: Decide what you want to do for each point + flagged by find_isolated_points.py, and edit + this script accordingly. Make sure the path to + the ROMS grid file is correct. Then open python + or ipython and type "run fix_isolated_points.py" + +romscice_sss_nudging.py: Interpolate the World Ocean Atlas 2013 monthly + climatology of sea surface salinity to the ROMS grid, + for use in surface salinity restoring. + To run: Make sure the paths to the ROMS grid file, + WOA input files, and desired output file + are correct (near the top of the file). Then + make sure that you have scipy version 0.14 or + higher (on raijin, this means switching to + python/2.7.6; instructions at the top of the + file). Then open python or ipython and type + "run romscice_sss_nudging.py". + ***POST-PROCESSING DIAGNOSTIC FILES*** @@ -287,6 +333,31 @@ common_grid.py: Interpolates ROMS output to a regular quarter-degree grid for CICE output files, and the desired output file on the common grid. +calc_rx1.py: Calculate the rx1 field (measure of grid steepness) for the given + grid and vertical stretching parameters. Save to a NetCDF file. + To run: Figure out what your vertical stretching parameters are + by looking at your *.in configuration file, or you can + just experiment based on other people's choices. Then + open python or ipython and type "run calc_rx1.py". The + script will prompt you for the paths to your ROMS grid + file and desired output file, as well as the vertical + stretching parameters. Note that this script assumes + Vtransform=2 (as do all scripts which call calc_z.py). + +plot_layers.py: Plot the terrain-following vertical levels through a given + longitude, in a latitude vs. depth slice. This is a good way + to test out different choices for vertical stretching + parameters. + To run: Figure out what your vertical stretching parameters are + by looking at your *.in configuration file, or you can + just experiment based on other people's choices. Then + open python or ipython and type "run plot_layers.py". + The script will prompt you for the ROMS grid file, the + longitude to interpolate to, and your vertical + stretching parameters. The plot will be displayed on + screen, and the script will repeat as many times as you + like. + ***DIAGNOSTIC FIGURES*** @@ -385,6 +456,29 @@ timeseries_3D.py: Calculates and plots timeseries of total ocean heat content, for the paths to the ROMS grid file, the ocean history or averages file, and the log file. +timeseries_sss.py: Calculates and plots timeseries of area-averaged sea surface + salinity, surface salt flux, and surface salt flux due to + salinity restoring during a ROMS simulation. Also writes the + timeseries to a log file so it doesn't have to be recomputed + if the run is extended; at the beginning of this script, + previous values will be read from the same log file if it + exists. + To run: Open python or ipython and type + "run timeseries_sss.py". The script will prompt you + for the path to the ROMS averages file and the log + file. + +timeseries_i2osalt.py: Calculates and plots timeseries of area-averaged sea ice + to ocean salt flux during a ROMS simulation. Also writes + the timeseries to a log file so it doesn't have to be + recomputed if the run is extended; at the beginning of + this script, previous values will be read from the same + log file if it exists. + To run: Open python or ipython and type + "run timeseries_i2osalt.py". The script will + prompt you for the path to the CICE history file + and the log file. + aice_animation.py: Create an animation of sea ice concentration for the given simulation, and save as an mp4 file. To run: If you are on raijin, first type "module load ffmpeg" @@ -943,27 +1037,18 @@ adv_u3_sss_anom.py: Plot the sea surface salinity anomaly between the U3 and The figure will be saved under the filename sss_u3_anom.png. -timeseries_sss.py: Calculates and plots timeseries of area-averaged sea surface - salinity and surface salt flux during a ROMS simulation. - Also writes the timeseries to a log file so it doesn't have - to be recomputed if the run is extended; at the beginning of - this script, previous values will be read from the same log - file if it exists. - To run: Open python or ipython and type - "run timeseries_sss.py". The script will prompt you - for the path to the ROMS averages file and the log - file. +adv_mld_2.py: Like adv_mld.py, but only showing U3_LIM and C4_LD. + To run: Make sure the paths to the simulation directories and 23 + August ROMS files are correct. Then open python or ipython + and type "run adv_mld_2.py". The figure will be saved + under the filename adv_mld_2.png. -timeseries_i2osalt.py: Calculates and plots timeseries of area-averaged sea ice - to ocean salt flux during a ROMS simulation. Also writes - the timeseries to a log file so it doesn't have to be - recomputed if the run is extended; at the beginning of - this script, previous values will be read from the same - log file if it exists. - To run: Open python or ipython and type - "run timeseries_i2osalt.py". The script will - prompt you for the path to the CICE history file - and the log file. +adv_thickness_2.py: Like adv_thickness.py, but only showing U3_LIM and C4_LD. + To run: Make sure the paths to the simulation directories + and 23 August ROMS files are correct. Then open + python or ipython and type "run adv_thickness_2.py". + The figure will be saved under the filename + adv_thickness_2.png. ***FIGURES FOR GARY'S LIMITERS PAPER*** diff --git a/find_isolated_points.py b/find_isolated_points.py index 74a82fb..5b60e63 100644 --- a/find_isolated_points.py +++ b/find_isolated_points.py @@ -1,28 +1,42 @@ from netCDF4 import Dataset from numpy import * +# Find all of the CICE grid points which are land (or ice shelf) on 3 sides. +# Sea ice can grow in these isolated points but cannot escape due to CICE's +# coastal boundary conditions, so it gets crazy thick (like 2 km thick). +# Print the indices of these points to the screen. This script assumes a +# periodic boundary in the longitude direction. +# Input: cice_kmt_file = path to CICE land mask file, created using cice_grid.py def find_isolated_points (cice_kmt_file): + # j-indices that might have sea ice start_j = 50 end_j = 250 + # Read land mask id = Dataset(cice_kmt_file, 'r') kmt = id.variables['kmt'][:,:] id.close() num_i = size(kmt,1) + # Double loop, can't find a cleaner way to do this for j in range(start_j, end_j): - for i in range(num_i-1): + for i in range(num_i): + # Check for unmasked points if kmt[j,i] == 1: + # Count the number of neighbours which are unmasked if i == num_i-1: + # Loop back to the beginning for neighbour on the left neighbours = array([kmt[j,i-1], kmt[j,0], kmt[j-1,i], kmt[j+1,i]]) else: neighbours = array([kmt[j,i-1], kmt[j,i+1], kmt[j-1,i], kmt[j+1,i]]) if sum(neighbours) < 2: + # Blocked on at least 3 sides print "i=" + str(i+1) + ', j=' + str(j+1) +# Command-line interface if __name__ == "__main__": cice_kmt_file = raw_input("Path to CICE land mask file: ") diff --git a/fix_isolated_points.py b/fix_isolated_points.py index efbdf72..5d3f5cc 100644 --- a/fix_isolated_points.py +++ b/fix_isolated_points.py @@ -1,10 +1,21 @@ from netCDF4 import Dataset from numpy import * +# Based on the results of find_isolated_points.py, fix all the points which +# are isolated on 3 sides in the CICE grid, by editing the ROMS grid. Set some +# of them to ice shelf points (with zice the average of the correct neighbour +# values), some to land points, remove the ice shelves on others, and make other +# land points ocean points (with h the average of the correct neighbour values). +# After running this script, rerun cice_grid.py to update the CICE grid. +# This script is grid-specific and manually written. If you are editing it +# for a new grid, make sure you take into account the ghost cells in ROMS +# and CICE: point (i,j) in CICE is point (i+1,j+1) in ROMS. def fix_isolated_pts (): + # Path to ROMS grid file grid_file = '../metroms_iceshelf/apps/common/grid/circ30S_quarterdegree_tmp.nc' + # Read the relevant fields id = Dataset(grid_file, 'a') mask_rho = id.variables['mask_rho'][:,:] h = id.variables['h'][:,:] # Fill value 50 @@ -165,6 +176,8 @@ def fix_isolated_pts (): id.variables['zice'][:,:] = zice id.close() + +# Command-line interface if __name__ == "__main__": fix_isolated_pts() diff --git a/plot_layers.py b/plot_layers.py index deec5ad..39c8f86 100644 --- a/plot_layers.py +++ b/plot_layers.py @@ -4,8 +4,18 @@ from calc_z import * from interp_lon_roms import * +# Plot the terrain-following vertical levels through a given longitude, in a +# latitude vs. depth slice. This is a good way to test out different choices +# for vertical stretching parameters. +# Input: +# grid_path = path to ROMS grid file +# lon0 = longitude to interpolate to (-180 to 180) +# depth_min = deepest depth to plot (negative, in metres) +# Vstretching, theta_s, theta_b, hc, N = vertical stretching parameters (see +# *.in configuration file if unsure) def plot_layers (grid_path, lon0, depth_min, Vstretching, theta_s, theta_b, hc, N): + # Read grid id = Dataset(grid_path, 'r') lon_2d = id.variables['lon_rho'][:,:] lat_2d = id.variables['lat_rho'][:,:] @@ -14,22 +24,29 @@ def plot_layers (grid_path, lon0, depth_min, Vstretching, theta_s, theta_b, hc, mask = id.variables['mask_rho'][:,:] id.close() + # Calculate a 3D field of depth values z_3d, sc_r, Cs_r = calc_z(h, zice, theta_s, theta_b, hc, N, None, Vstretching) + # Figure out how to write longitude in the title if lon0 < 0: lon_string = str(int(round(-lon0))) + r'$^{\circ}$W' else: lon_string = str(int(round(lon0))) + r'$^{\circ}$E' + # Get longitude in the range (0,360) as per ROMS convention if lon0 < 0: lon0 += 360 + # Get a 3D land mask mask_3d = tile(mask, (N,1,1)) + # Interpolate land mask, depth, longitude to the given longitude mask_interp, z, lat = interp_lon_roms(mask_3d, z_3d, lat_2d, lon_2d, lon0) + # Simple array containing layer numbers, same size as z layer = zeros(shape(z)) for k in range(N): layer[k,:] = k+1 + # Mask out land layer = ma.masked_where(mask_interp==0, layer) # Choose latitude bounds based on land mask @@ -58,8 +75,10 @@ def plot_layers (grid_path, lon0, depth_min, Vstretching, theta_s, theta_b, hc, #lat_min = -85 #lat_max = -75 + # Contour levels lev = range(1, N) + # Plot fig = figure(figsize=(18,6)) #(18,12)) contour(lat, z, layer, lev, colors='k') title(lon_string) #(r"ROMS terrain-following levels through the Ross Ice Shelf cavity (180$^{\circ}$E)", fontsize=24) @@ -72,10 +91,12 @@ def plot_layers (grid_path, lon0, depth_min, Vstretching, theta_s, theta_b, hc, #fig.savefig('ross_levels.png') fig.show() + # Put longitude back to original range in case the user wants to go again if lon0 > 180: lon0 -= 360 +# Command-line interface if __name__ == "__main__": grid_path = raw_input("Path to grid file: ") @@ -88,6 +109,7 @@ def plot_layers (grid_path, lon0, depth_min, Vstretching, theta_s, theta_b, hc, N = int(raw_input("N: ")) plot_layers (grid_path, lon0, depth_min, Vstretching, theta_s, theta_b, hc, N) + # Keep repeating until the user wants to exit while True: repeat = raw_input("Make another plot (y/n)? ") if repeat == 'y': diff --git a/romscice_sss_nudging.py b/romscice_sss_nudging.py index 79e10ac..10a4fcd 100644 --- a/romscice_sss_nudging.py +++ b/romscice_sss_nudging.py @@ -3,20 +3,36 @@ from scipy.interpolate import RegularGridInterpolator from scipy.spatial import KDTree +# Interpolate the World Ocean Atlas 2013 monthly climatology of sea surface +# salinity to the ROMS grid, for use in surface salinity restoring. + +# 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 + def romscice_sss_nudging (): + # Path to ROMS grid file grid_file = '../metroms_iceshelf/apps/common/grid/circ30S_quarterdegree.nc' + # Paths to reconstruct WOA 2013 monthly files woa_head = '../metroms_iceshelf/data/originals/WOA_2013/woa13_95A4_s' woa_tail = '_01v2.nc' + # Path to output file output_file = '../metroms_iceshelf/data/sss_nudge.nc' + # Northernmost j-index of WOA grid to read nbdry_woa = 61 - print 'Reading ECCO2 grid' + # Read WOA grid + print 'Reading WOA grid' id = Dataset(woa_head + '01' + woa_tail, 'r') lon_woa_raw = id.variables['lon'][:] lat_woa = id.variables['lat'][0:nbdry_woa] id.close() + # Make sure longitude goes from -180 to 180 index = lon_woa_raw > 180 lon_woa_raw[index] -= 360 @@ -41,9 +57,11 @@ def romscice_sss_nudging (): num_lon = size(lon_roms, 1) num_lat = size(lon_roms, 0) + # Make sure longitude goes from -180 to 180 index = lon_roms > 180 lon_roms[index] -= 360 + # Set up output file print 'Setting up ' + output_file out_id = Dataset(output_file, 'w') out_id.createDimension('xi_rho', num_lon) @@ -57,35 +75,57 @@ def romscice_sss_nudging (): out_id.variables['SSS'].long_name = 'surface salinity' out_id.variables['SSS'].units = 'psu' + # Loop over months for month in range(12): print 'Processing month ' + str(month+1) + ' of 12' + # Construct file path if month+1 < 10: filename = woa_head + '0' + str(month+1) + woa_tail else: filename = woa_head + str(month+1) + woa_tail + # Read surface salinity id = Dataset(filename, 'r') sss_raw = transpose(id.variables['s_an'][0,0,0:nbdry_woa,:]) id.close() + # Wrap the longitude axis sss = ma.array(zeros((size(lon_woa), size(lat_woa)))) sss[1:-1,:] = ma.copy(sss_raw) sss[0,:] = ma.copy(sss_raw[-1,:]) sss[-1,:] = ma.copy(sss_raw[0,:]) + # Interpolate to ROMS grid sss_interp = interp_woa2roms_sfc(sss, lon_woa, lat_woa, lon_roms, lat_roms) + # Calculate time value: 12 values equally spaced throughout the year time = 365.25/12*(month+0.5) + # Save to file out_id.variables['sss_time'][month] = time out_id.variables['SSS'][month,:,:] = sss_interp out_id.close() +# Given a surface variable on the World Ocean Atlas grid, fill the land mask +# with nearest neighbours and interpolate onto the ROMS grid. +# Input: +# A = array of size nxmxo containing values on the World Ocean Atlas grid +# (dimension longitude x latitude) +# lon_woa = array of length n containing WOA longitude values +# lat_woa = array of length m containing WOA latitude values +# lon_roms = array of size pxq containing ROMS longitude values +# lat_roms = array of size pxq containing ROMS latitude values +# Output: +# B = array of size pxq containing values on the ROMS grid (dimension latitude +# x longitude) def interp_woa2roms_sfc (A, lon_woa, lat_woa, lon_roms, lat_roms): + # Fill land mask with nearest neighbours j,i = mgrid[0:A.shape[0], 0:A.shape[1]] jigood = array((j[~A.mask], i[~A.mask])).T jibad = array((j[A.mask], i[A.mask])).T A[A.mask] = A[~A.mask][KDTree(jigood).query(jibad)[1]] + # Build a function for linear interpolation of A interp_function = RegularGridInterpolator((lon_woa, lat_woa), A) + # Call this function for each point on the ROMS grid (vectorised) B = interp_function((lon_roms, lat_roms)) return B