Bug fixes and some new scripts for diagnosing sea ice weirdness

knaughten · Jan 23, 2017 · f3ad52c · f3ad52c
1 parent 6c2a760
commit f3ad52c
Show file tree

Hide file tree

Showing 20 changed files with 1,348 additions and 103 deletions.
diff --git a/adv_amery_tsplots.py b/adv_amery_tsplots.py
@@ -3,56 +3,80 @@
 from matplotlib.pyplot import *
 from calc_z import *
 
+# Create a 2x2 plot showing zonal slices of temperature and salinity through
+# 71E (Amery Ice Shelf) at the end of the U3_LIM and C4_LD simulations.
 def adv_amery_tsplots ():
 
+    # Paths to simulation directories
     paths = ['/short/m68/kaa561/advection/u3_lim/', '/short/m68/kaa561/advection/c4_l/']
+    # Titles for plotting
     labels = [r'a) Temperature ($^{\circ}$C), U3_LIM', r'b) Temperature ($^{\circ}$C), C4_LD', 'c) Salinity (psu), U3_LIM', 'd) Salinity (psu), C4_LD']
+    # File name: daily average for 31 December
     file_tail = 'ocean_avg_31dec.nc'
     var_names = ['temp', 'salt']
-    tstep = 1 #366
+    # If 31 December doesn't have its own file, put the time index here
+    tstep = 1 #366 if all one file of daily averages for entire simulation
+    # Longitude to interpolate to
     lon0 = 71
+    # Deepest depth to plot
     depth_min = -500
+    # Bounds and ticks for colour scales
     scale_min = [-2, 33.8]
-    scale_max = [3, 34.8]
+    scale_max = [3, 34.8]    
     scale_ticks = [1, 0.2]
+    # Bounds on latitude
     lat_min = -72
     lat_max = -50
-
+    # Grid parameters
     theta_s = 4.0
     theta_b = 0.9
     hc = 40
     N = 31
 
+    # Set up figure
     fig = figure(figsize=(18,12))
+    # Loop over simulations
     for sim in range(2):
+        # Loop over variables (temp and salt)
         for var in range(2):
+            # Read 3D variable
             id = Dataset(paths[sim]+file_tail, 'r')
             data_3d = id.variables[var_names[var]][tstep-1,:,:,:]
+            # Also read sea surface height
             zeta = id.variables['zeta'][tstep-1,:,:]
             if sim==0 and var==0:
+                # For the first simulation, read grid variables
                 h = id.variables['h'][:,:]
                 zice = id.variables['zice'][:,:]
                 lon_2d = id.variables['lon_rho'][:,:]
                 lat_2d = id.variables['lat_rho'][:,:]
             id.close()
+            # Calculate 3D z-coordinates
             z_3d, sc_r, Cs_r = calc_z(h, zice, theta_s, theta_b, hc, N, zeta)
+            # Interpolate temp/salt, z-coordinates, and latitude to 71E
             data, z, lat = interp_lon(data_3d, z_3d, lat_2d, lon_2d, lon0)
             ax = fig.add_subplot(2, 2, 2*var+sim+1)
+            # Shade data (pcolor not contourf so we don't misrepresent the
+            # model grid)
             img = pcolor(lat, z, data, vmin=scale_min[var], vmax=scale_max[var], cmap='jet')
+            # Add title
             title(labels[2*var+sim], fontsize=24)
+            # Label axes
             if var == 1:
                 xlabel('Latitude', fontsize=16)
             if sim == 0:
                 ylabel('Depth (m)', fontsize=16)
             xlim([lat_min, lat_max])
             ylim([depth_min, 0])
             if sim == 1:
+                # Add colorbars for each variable
                 if var == 0:
                     cbaxes = fig.add_axes([0.93, 0.575, 0.01, 0.3])
                 elif var == 1:
                     cbaxes = fig.add_axes([0.93, 0.125, 0.01, 0.3])
                 cbar = colorbar(img, ticks=arange(scale_min[var], scale_max[var]+scale_ticks[var], scale_ticks[var]), cax=cbaxes, extend='both')
                 cbar.ax.tick_params(labelsize=14)
+            # Set ticks the way we want them
             lat_ticks = arange(lat_min+2, lat_max+1, 5)
             ax.set_xticks(lat_ticks)
             lat_labels = []
@@ -66,6 +90,7 @@ def adv_amery_tsplots ():
                 depth_labels.append(str(int(round(-val))))
             ax.set_yticklabels(depth_labels, fontsize=14)
 
+    # Main title
     suptitle(r'71$^{\circ}$E (Amery Ice Shelf), 31 December', fontsize=30)
     #fig.show()
     fig.savefig('adv_amery_tsplots.png')
@@ -176,6 +201,7 @@ def interp_lon_helper (lon, lon0):
     return ie, iw, coeff1, coeff2
 
 
+# Command-line interface
 if __name__ == "__main__":
 
     adv_amery_tsplots()
diff --git a/adv_frazil.py b/adv_frazil.py
@@ -3,30 +3,43 @@
 from matplotlib.pyplot import *
 #import colormaps as cmaps
 
+# Create a 3x2 plot of annually averaged frazil formation for each advection
+# experiment: the absolute values for U3_LIM, and the anomalies from U3_LIM for
+# the other 5 experiments.
 def adv_frazil ():
 
+    # Paths to simulation directories
     paths = ['/short/m68/kaa561/advection/u3_lim/', '/short/m68/kaa561/advection/u3/', '/short/m68/kaa561/advection/c4_l/', '/short/m68/kaa561/advection/c4_h/', '/short/m68/kaa561/advection/a4_l/', '/short/m68/kaa561/advection/a4_h/']
+    # Titles for plotting
     labels = ['a) U3_LIM', 'b) U3 - U3_LIM', 'c) C4_LD - U3_LIM', 'd) C4_HD - U3_LIM', 'e) A4_LD - U3_LIM', 'f) A4_HD - U3_LIM']
+    # File name: annual average
     file_tail = 'iceh_avg.nc'
 
+    # Bounds and ticks for colour scales
     max_abs = 0.2
     tick_abs = 0.05
     max_anom = 1.0
     tick_anom = 0.5
 
+    # Degrees to radians conversion factor
     deg2rad = pi/180.
+    # Centre of missing circle in grid
     lon_c = 50
     lat_c = -83
+    # Radius of missing circle
     radius = 10.5
+    # Boundary of regular grid to embed circle in
     circle_bdry = -70+90
-
+
+    # Read frazil data from U3_LIM simulation; also grid and mask variables
     id = Dataset(paths[0]+file_tail, 'r')
     data_tmp = id.variables['frazil'][0,:350,:]
     lon_tmp = id.variables['TLON'][:350,:]
     lat_tmp = id.variables['TLAT'][:350,:]
     mask_tmp = id.variables['tmask'][:350,:]
     id.close()
 
+    # Wrap periodic boundary so there isn't a gap in the plot
     lon = ma.empty([size(lon_tmp,0), size(lon_tmp,1)+1])
     lat = ma.empty([size(lat_tmp,0), size(lat_tmp,1)+1])
     mask = ma.empty([size(mask_tmp,0), size(mask_tmp,1)+1])
@@ -40,53 +53,73 @@ def adv_frazil ():
     data0[:,:-1] = data_tmp
     data0[:,-1] = data_tmp[:,0]
 
+    # Land mask
     land = ma.masked_where(mask==1, mask)
 
+    # Circumpolar x and y coordinates for plotting
     x = -(lat+90)*cos(lon*deg2rad+pi/2)
     y = (lat+90)*sin(lon*deg2rad+pi/2)
+    # Coordinates of centre of missing circle
     x_c = -(lat_c+90)*cos(lon_c*deg2rad+pi/2)
     y_c = (lat_c+90)*sin(lon_c*deg2rad+pi/2)
+    # Regular grid to embed missing circle in
     x_reg, y_reg = meshgrid(linspace(-circle_bdry, circle_bdry, num=100), linspace(-circle_bdry, circle_bdry, num=100))
+    # Mask everything except the circle out of the regular grid
     land_circle = zeros(shape(x_reg))
     land_circle = ma.masked_where(sqrt((x_reg-x_c)**2 + (y_reg-y_c)**2) > radius, land_circle)
 
+    # Set up figure
     fig = figure(figsize=(9,15))
     ax = fig.add_subplot(3, 2, 1, aspect='equal')
+    # First shade land
     contourf(x, y, land, 1, colors=(('0.6', '0.6', '0.6')))
+    # Fill in missing circle
     contourf(x_reg, y_reg, land_circle, 1, colors=(('0.6', '0.6', '0.6')))
+    # Shade the frazil data (pcolor not contourf so we don't misrepresent the
+    # model grid)
     img0 = pcolor(x, y, data0, vmin=0, vmax=max_abs, cmap='jet')
     axis('off')
+    # Add title
     title(labels[0], fontsize=20)
+    # Add colorbar
     cbaxes0 = fig.add_axes([0.025, 0.7, 0.02, 0.2])
     cbar0 = colorbar(img0, ticks=arange(0, max_abs+tick_abs, tick_abs), cax=cbaxes0, extend='max')
     cbar0.ax.tick_params(labelsize=16)  
 
+    # Loop over the other simulations
     for sim in range(1, len(paths)):
+        # Read the frazil data
         id = Dataset(paths[sim]+file_tail, 'r')
         data_tmp = id.variables['frazil'][0,:350,:]
         id.close()
-        data = ma.empty([size(data_tmp,0), size(data_tmp,1)+1])
+        # Wrap the periodic boundary
+        data = ma.empty([size(data_tmp,0), size(data_tmp,1)+1])        
         data[:,:-1] = data_tmp
         data[:,-1] = data_tmp[:,0]
+        # Calculate anomaly from U3_LIM
         data = data - data0
+        # Add to plot, same as before
         ax = fig.add_subplot(3, 2, sim+1, aspect='equal')
         contourf(x, y, land, 1, colors=(('0.6', '0.6', '0.6')))
         contourf(x_reg, y_reg, land_circle, 1, colors=(('0.6', '0.6', '0.6')))
         img = pcolor(x, y, data, vmin=-max_anom, vmax=max_anom, cmap='RdBu_r')
         axis('off')
         title(labels[sim], fontsize=20)
         if sim == 3:
+            # Only add an anomaly colorbar for one of the simulations
             cbaxes = fig.add_axes([0.025, 0.4, 0.02, 0.2])
             cbar = colorbar(img, ticks=arange(-max_anom, max_anom+tick_anom, tick_anom), cax=cbaxes, extend='both')
             cbar.ax.tick_params(labelsize=16)
-
+
+    # Main title
     suptitle('Annually averaged frazil ice formation (cm/day)', fontsize=24)
     subplots_adjust(wspace=0.025,hspace=0.2)
 
     #fig.show()
     fig.savefig('adv_frazil.png')
 
 
+# Command-line interface
 if __name__ == '__main__':
 
     adv_frazil()

diff --git a/adv_mld.py b/adv_mld.py
@@ -3,97 +3,116 @@
 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.
 def adv_mld ():
 
+    # Paths to simulation directories
     paths = ['/short/m68/kaa561/advection/u3_lim/', '/short/m68/kaa561/advection/u3/', '/short/m68/kaa561/advection/c4_l/', '/short/m68/kaa561/advection/c4_h/', '/short/m68/kaa561/advection/a4_l/', '/short/m68/kaa561/advection/a4_h/']
+    # Titles for plotting
     labels = ['a) U3_LIM', 'b) U3 - U3_LIM', 'c) C4_LD - U3_LIM', 'd) C4_HD - U3_LIM', 'e) A4_LD - U3_LIM', 'f) A4_HD - U3_LIM']
+    # File name: daily average for 23 August
     file_tail = 'ocean_avg_23aug.nc'
-    tstep = 1 #236
+    # If 23 August doesn't have its own file, put the time index here
+    tstep = 1 #236 if all one file of daily averages for entire simulation
 
+    # Bounds and ticks for colour scales
     max_abs = 300
     tick_abs = 100
     max_anom = 200
     tick_anom = 100
 
+    # Degrees to radians conversion factor
     deg2rad = pi/180.
+    # Centre of missing circle in grid
     lon_c = 50
     lat_c = -83
+    # Radius of missing circle
     radius = 10.5
+    # Boundary of regular grid to embed circle in
     circle_bdry = -70+90
-
+
+    # Read mixed layer depth from U3_LIM simulation; also grid and mask
+    # variables
     id = Dataset(paths[0]+file_tail, 'r')
-    data_tmp = id.variables['Hsbl'][tstep-1,:350,:]
-    lon_tmp = id.variables['lon_rho'][:350,:]
-    lat_tmp = id.variables['lat_rho'][:350,:]
-    mask_tmp = id.variables['mask_rho'][:350,:]
-    zice_tmp = id.variables['zice'][:350,:]
+    data = id.variables['Hsbl'][tstep-1,:350,1:]
+    lon = id.variables['lon_rho'][:350,1:]
+    lat = id.variables['lat_rho'][:350,1:]
+    mask = id.variables['mask_rho'][:350,1:]
+    zice = id.variables['zice'][:350,1:]
     id.close()
 
-    index = zice_tmp != 0
-    mask_tmp[index] = 0.0
-    data_tmp = ma.masked_where(zice_tmp!=0, -data_tmp)
-
-    lon = ma.empty([size(lon_tmp,0), size(lon_tmp,1)+1])
-    lat = ma.empty([size(lat_tmp,0), size(lat_tmp,1)+1])
-    mask = ma.empty([size(mask_tmp,0), size(mask_tmp,1)+1])
-    data0 = ma.empty([size(data_tmp,0), size(data_tmp,1)+1])
-    lon[:,:-1] = lon_tmp
-    lon[:,-1] = lon_tmp[:,0]
-    lat[:,:-1] = lat_tmp
-    lat[:,-1] = lat_tmp[:,0]
-    mask[:,:-1] = mask_tmp
-    mask[:,-1] = mask_tmp[:,0]
-    data0[:,:-1] = data_tmp
-    data0[:,-1] = data_tmp[:,0]
+    # Mask out the ice shelf cavities and switch sign on mixed layer depth
+    index = zice != 0
+    mask[index] = 0.0
+    data = ma.masked_where(zice!=0, -data)
 
+    # Land mask
     land = ma.masked_where(mask==1, mask)
 
+    # Circumpolar x and y coordinates for plotting
     x = -(lat+90)*cos(lon*deg2rad+pi/2)
     y = (lat+90)*sin(lon*deg2rad+pi/2)
+    # Coordinates of centre of missing circle
     x_c = -(lat_c+90)*cos(lon_c*deg2rad+pi/2)
     y_c = (lat_c+90)*sin(lon_c*deg2rad+pi/2)
+    # Regular grid to embed missing circle in 
     x_reg, y_reg = meshgrid(linspace(-circle_bdry, circle_bdry, num=100), linspace(-circle_bdry, circle_bdry, num=100))
+    # Mask everything except the circle out of the regular grid
     land_circle = zeros(shape(x_reg))
     land_circle = ma.masked_where(sqrt((x_reg-x_c)**2 + (y_reg-y_c)**2) > radius, land_circle)
 
+    # Set up figure
     fig = figure(figsize=(9,15))
     ax = fig.add_subplot(3, 2, 1, aspect='equal')
+    # First shade land
     contourf(x, y, land, 1, colors=(('0.6', '0.6', '0.6')))
+    # Fill in missing circle
     contourf(x_reg, y_reg, land_circle, 1, colors=(('0.6', '0.6', '0.6')))
+    # Shade the mixed layer depth (pcolor not contourf so we don't misrepresent
+    # the model grid)
     img0 = pcolor(x, y, data0, vmin=0, vmax=max_abs, cmap='jet') #cmaps.viridis)
     axis('off')
+    # Add title
     title(labels[0], fontsize=20)
+    # Add colorbar
     cbaxes0 = fig.add_axes([0.025, 0.7, 0.02, 0.2])
     cbar0 = colorbar(img0, ticks=arange(0, max_abs+tick_abs, tick_abs), cax=cbaxes0, extend='max')
     cbar0.ax.tick_params(labelsize=16)  
 
+    # Loop over the other simulations
     for sim in range(1, len(paths)):
+        # Read mixed layer depth
         id = Dataset(paths[sim]+file_tail, 'r')
-        data_tmp = id.variables['Hsbl'][tstep-1,:350,:]
+        data = id.variables['Hsbl'][tstep-1,:350,:]
         id.close()
-        data_tmp = ma.masked_where(zice_tmp!=0, -data_tmp)
-        data = ma.empty([size(data_tmp,0), size(data_tmp,1)+1])
-        data[:,:-1] = data_tmp
-        data[:,-1] = data_tmp[:,0]
+        # Mask out the ice shelf cavities and switch sign
+        data = ma.masked_where(zice!=0, -data)
+        # Calculate anomaly from U3_LIM
         data = data - data0
+        # Add to plot, same as before
         ax = fig.add_subplot(3, 2, sim+1, aspect='equal')
         contourf(x, y, land, 1, colors=(('0.6', '0.6', '0.6')))
         contourf(x_reg, y_reg, land_circle, 1, colors=(('0.6', '0.6', '0.6')))
         img = pcolor(x, y, data, vmin=-max_anom, vmax=max_anom, cmap='RdBu_r')
         axis('off')
         title(labels[sim], fontsize=20)
         if sim == 3:
+            # Only add an anomaly colorbar for one of the simulations
             cbaxes = fig.add_axes([0.025, 0.4, 0.02, 0.2])
             cbar = colorbar(img, ticks=arange(-max_anom, max_anom+tick_anom, tick_anom), cax=cbaxes, extend='both')
             cbar.ax.tick_params(labelsize=16)
-
+
+    # Main title
     suptitle('Mixed layer depth (m) on 23 August', fontsize=28)
     subplots_adjust(wspace=0.025,hspace=0.15)
 
     #fig.show()
     fig.savefig('adv_mld.png')
 
 
+# Command-line interface
 if __name__ == '__main__':
 
     adv_mld()