Comments and documentation

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

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

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,31 +19,40 @@ 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))
     id.createVariable('rx1', 'f8', ('eta_rho', 'xi_rho'))
     id.variables['rx1'][:,:] = rx1
     id.close()
 
-
+
+# Command-line interface
 if __name__ == "__main__":
 
     grid_path = raw_input("Path to grid file: ")

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***

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: ")

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

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':