Documentation! Finally!

add_tide_period.py

@@ -8,7 +8,7 @@ def add_tide_period ():
 
     # Tide file to add periods to
     # Order is q1, o1, p1, k1, n2, m2, s2, k2
-    tide_file = '../ROMS-CICE-MCT/data/pot_tides.nc'
+    tide_file = '../metroms_iceshelf/data/pot_tides.nc'
     # Periods in seconds
     # Order is m2, s2, n2, k2, k1, o1, p1, q1, mf, mm
     period = array([44714.165191868, 43200.0012869521, 86164.0770050671, 92949.6357005365, 45570.0535117177, 86637.1997716528, 43082.0503185947, 96726.0857029666, 2380715.86358729, 1180295.54554976])

adv_polynyas.py

@@ -2,7 +2,7 @@
 from numpy import *
 from matplotlib.pyplot import *
 
-# Create a 1x2 plot of sea ice concentration along the coast of Wilkes Land
+# Create a 1x2 plot of sea ice concentration near the Amery Ice Shelf
 # on 23 August (the sea ice area max), showing the difference in polynya size
 # between the U3_LIM and C4_LD simulations.
 def adv_polynyas ():

adv_ross_tsplots.py

@@ -93,8 +93,8 @@ def adv_ross_tsplots ():
 
     # Main title
     suptitle(r'180$^{\circ}$E, 31 December (Ross Sea)', fontsize=30)
-    fig.show()
-    #fig.savefig('adv_ross_tsplots.png')
+    #fig.show()
+    fig.savefig('adv_ross_tsplots.png')
 
 
 # Command-line interface

salinity_diff.py → adv_u3_sss_anom.py

@@ -2,10 +2,15 @@
 from numpy import *
 from matplotlib.pyplot import *
 
-def salinity_diff ():
+# Plot the sea surface salinity anomaly between the U3 and U3_LIM simulations
+# on 23 August (sea ice area maximum).
+def adv_u3_sss_anom ():
 
+    # Paths to simulation directories
     paths = ['/short/m68/kaa561/advection/u3_lim/', '/short/m68/kaa561/advection/u3/']
+    # File for 23 August daily average
     file_tail = 'iceh.1992-08-23.nc'
+    # Bounds on colour scale
     max_anom = 0.2
     tick_anom = 0.1
     # Degrees to radians conversion factor
@@ -75,8 +80,9 @@ def salinity_diff ():
     fig.savefig('sss_u3_anom.png')
 
 
+# Command-line interface
 if __name__ == "__main__":
 
-    salinity_diff()
+    adv_u3_sss_anom()
 
 

cice_ini.py

@@ -2,48 +2,72 @@
 from numpy import *
 from scipy.interpolate import griddata
 
+# Make a CICE initial conditions file (non-standard option ice_ic='mask') based
+# on NSIDC satellite observations of sea ice concentration. Wherever 
+# concentration exceeds 15%, set it to 100%; everywhere else, set it to 0%. 
+# The CICE code will read this "mask" and set a constant sea ice thickness of 
+# 1 m and snow thickness of 0.2 m.
+# Input:
+# nsidc_file = path to NSIDC file: CDR monthly average of sea ice concentration
+# roms_grid = path to ROMS grid file
+# out_file = desired path to output file
 def cice_ini (nsidc_file, roms_grid, out_file):
 
+    # Read the NSIDC grid and data
     id = Dataset(nsidc_file, 'r')
     nsidc_lon = id.variables['longitude'][:,:]
     nsidc_lat = id.variables['latitude'][:,:]
     nsidc_aice = id.variables['seaice_conc_monthly_cdr'][0,:,:]
     id.close()
 
+    # Make sure longitude goes from 0 to 360 to fit with ROMS convention
     index = nsidc_lon < 0
     nsidc_lon[index] = nsidc_lon[index] + 360
 
+    # Make sure sea ice concentration values go from 0 to 1
     index = nsidc_aice < 0
     nsidc_aice[index] = 0.0
     index = nsidc_aice > 1
     nsidc_aice[index] = 1.0
 
+    # Read the ROMS grid
     id = Dataset(roms_grid, 'r')
     roms_lon = id.variables['lon_rho'][:,:]
     roms_lat = id.variables['lat_rho'][:,:]
     id.close()
 
+    # Make an nx2 array of the NSIDC coordinates, flattened
     points = empty([size(nsidc_lon), 2])
     points[:,0] = ravel(nsidc_lon)
     points[:,1] = ravel(nsidc_lat)
+    # Also flatten the NSIDC data
     values = ravel(nsidc_aice)
+    # Now make an mx2 array of the ROMS coordinates
     xi = empty([size(roms_lon), 2])
     xi[:,0] = ravel(roms_lon)
     xi[:,1] = ravel(roms_lat)
+    # Linear interpolation to the ROMS coordinates
     result = griddata(points, values, xi, method='linear', fill_value=-999)
+    # Also do a nearest-neighbour interpolation
     result2 = griddata(points, values, xi, method='nearest')
+    # Replace missing values (from land mask) in linear interpolation with
+    # nearest-neighbour interpolation so there are no coastal artifacts
     index = result==-999
     result[index] = result2[index]
+    # Reshape to the correct dimensions
     aice = reshape(result, shape(roms_lon))
 
+    # Cut-off concentration of 15%
     index = aice > 0.15
     aice[index] = 1.0
     aice[~index] = 0.0
 
+    # Remove the edges (CICE only wants interior points)
     aice = aice[1:-1,1:-1]
     num_lon = size(aice, 1)
     num_lat = size(aice, 0)
 
+    # Write output
     out_id = Dataset(out_file, 'w')
     out_id.createDimension('xi_rho', num_lon)
     out_id.createDimension('eta_rho', num_lat)
@@ -58,6 +82,7 @@ def cice_ini (nsidc_file, roms_grid, out_file):
     out_id.close()
 
 
+# Command-line interface
 if __name__ == "__main__":
 
     nsidc_file = '/short/m68/kaa561/nsidc_aice/seaice_conc_monthly_sh_f11_199201_v02r00.nc'