Fixed vector rotation issue, fixed some grid issues, and a bunch of n…

…ew figure scripts.

add_iceberg_melt.py

@@ -12,10 +12,10 @@
 def add_iceberg_melt (file):
 
     # Naming conventions for iceberg files
-    iceberg_head = '/short/m68/kaa561/ROMS-CICE-MCT/data/MartinAdcroft2010_iceberg_meltfluxes/icebergs.1861-1960.'
+    iceberg_head = '/short/m68/kaa561/ROMS-CICE-MCT/data/originals/MartinAdcroft2010_iceberg_meltfluxes/icebergs.1861-1960.'
     iceberg_tail = '.melt.nc'
     # Iceberg grid file
-    iceberg_grid = '/short/m68/kaa561/ROMS-CICE-MCT/data/MartinAdcroft2010_iceberg_meltfluxes/icebergs.static.nc'
+    iceberg_grid = '/short/m68/kaa561/ROMS-CICE-MCT/data/originals/MartinAdcroft2010_iceberg_meltfluxes/icebergs.static.nc'
     # ROMS grid file
     roms_grid ='/short/m68/kaa561/ROMS-CICE-MCT/apps/common/grid/circ30S_quarterdegree_10m.nc'
     # Density of freshwater
@@ -58,7 +58,7 @@ def add_iceberg_melt (file):
         melt_roms = melt_roms/rho_w*seconds_per_12h
         # Add to precipitation field for this month
         id = Dataset(file, 'a')
-        id.variables['rain'][month,:,:] = id.variables['rain'][month,:,:] - melt_roms
+        id.variables['rain'][month,:,:] = id.variables['rain'][month,:,:] + melt_roms
         id.close()
 
 

add_tide_period.py

@@ -0,0 +1,30 @@
+from netCDF4 import Dataset
+from numpy import *
+
+# Given a ROMS tide file with tidal potential amplitude and phase angle
+# (created using Kate's potential tide scripts), add the tidal period in seconds
+# for each component.
+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'
+    # 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])
+    # Put them in the right order
+    period_reorder = array([period[7], period[5], period[6], period[4], period[2], period[0], period[1], period[3]])
+
+    id = Dataset(tide_file, 'a')
+    id.createVariable('tide_period', 'f8', ('tide_period'))
+    id.variables['tide_period'].long_name = 'tide angular period'
+    id.variables['tide_period'].units = 'seconds'
+    id.variables['tide_period'][:] = period_reorder
+    id.close()
+
+
+# Command-line interface
+if __name__ == "__main__":
+
+    add_tide_period()
+

aice_animation.py

@@ -13,9 +13,9 @@
 # Directory containing CICE output files
 directory = '/short/y99/kaa561/roms_spinup_newest/cice/'
 # Number of time indices in each file
-num_ts = [504]
+num_ts = [270, 252, 252, 252, 252]
 # File number to start with for the animation (1-based)
-start_file = 1
+start_file = 5
 # Degrees to radians conversion factor
 deg2rad = pi/180
 # Names of each month for making titles
@@ -65,6 +65,6 @@ def animate(i):
     return img
 
 # Animate once every time index from start_file to the last file
-anim = FuncAnimation(fig, func=animate, frames=range(431,504)) #sum(array(num_ts[start_file-1:])))
+anim = FuncAnimation(fig, func=animate, frames=range(179,252)) #sum(array(num_ts[start_file-1:])))
 # Save as an mp4 with one frame per second
 anim.save('aice.mp4', fps=1)

bw_diff_plot.py

@@ -0,0 +1,99 @@
+from netCDF4 import Dataset
+from numpy import *
+from matplotlib.pyplot import *
+
+def bw_diff_plot (file_path, var_name, colour_bounds=None, save=False, fig_name=None):
+
+    init_path = '../ROMS-CICE-MCT/data/woa_ini.nc'
+    deg2rad = pi/180
+    lon_c = 50
+    lat_c = -83
+    radius = 10.1
+
+    id = Dataset(init_path, 'r')
+    data1 = id.variables[var_name][0,0,:-15,:-2]
+    id.close()
+
+    id = Dataset(file_path, 'r')
+    lon = id.variables['lon_rho'][:-15,:-2]
+    lat = id.variables['lat_rho'][:-15,:-2]
+    data2 = id.variables[var_name][-1,0,:-15,:-2]
+    id.close()
+
+    data_diff = data2 - data1
+
+    # Convert grid to spherical coordinates
+    x = -(lat+90)*cos(lon*deg2rad+pi/2)
+    y = (lat+90)*sin(lon*deg2rad+pi/2)
+    # Find centre in spherical coordinates
+    x_c = -(lat_c+90)*cos(lon_c*deg2rad+pi/2)
+    y_c = (lat_c+90)*sin(lon_c*deg2rad+pi/2)
+    # Build a regular x-y grid and select the missing circle
+    x_reg, y_reg = meshgrid(linspace(amin(x), amax(x), num=1000), linspace(amin(y), amax(y), num=1000))
+    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)
+
+    if colour_bounds is not None:
+        # User has set bounds on colour scale
+        lev = linspace(colour_bounds[0], colour_bounds[1], num=50)
+        if colour_bounds[0] == -colour_bounds[1]:
+            # Bounds are centered on zero, so choose a blue-to-red colourmap
+            # centered on yellow
+            colour_map = 'RdYlBu_r'
+        else:
+            colour_map = 'jet'
+    else:
+        max_val = amax(abs(data_diff))
+        lev = linspace(-max_val, max_val, num=50)
+        colour_map = 'RdYlBu_r'
+
+    # Plot
+    fig = figure(figsize=(16,12))
+    ax = fig.add_subplot(1,1,1,aspect='equal')
+    fig.patch.set_facecolor('white')
+    # First shade everything in grey
+    contourf(x, y, ones(shape(data_diff)), 1, colors=(('0.6', '0.6', '0.6')))
+    # Fill in the missing circle
+    contourf(x_reg, y_reg, land_circle, 1, colors=(('0.6', '0.6', '0.6')))
+    # Now shade the temperature (land mask will remain grey)
+    contourf(x, y, data_diff, lev, cmap=colour_map, extend='both')
+    cbar = colorbar()
+    cbar.ax.tick_params(labelsize=20)
+    if var_name == 'temp':
+        title(r'Difference in bottom water temperature since initialisation ($^{\circ}$C)', fontsize=30)
+    elif var_name == 'salt':
+        title(r'Difference in bottom water salinity since initialisation (psu)', fontsize=30)
+    axis('off')
+
+    # Finished
+    if save:
+        fig.savefig(fig_name)
+    else:
+        fig.show()
+
+
+if __name__ == "__main__":
+
+    file_path = raw_input("Path to ocean history/averages file: ")
+    tmp = raw_input("Temperature (t) or salinity (s)? ")
+    if tmp == 't':
+        var_name = 'temp'
+    elif tmp == 's':
+        var_name = 'salt'
+    colour_bounds = None
+    get_bounds = raw_input("Set bounds on colour scale (y/n)? ")
+    if get_bounds == 'y':
+        lower_bound = float(raw_input("Lower bound: "))
+        upper_bound = float(raw_input("Upper bound: "))
+        colour_bounds = [lower_bound, upper_bound]
+    action = raw_input("Save figure (s) or display on screen (d)? ")
+    if action == 's':
+        save = True
+        fig_name = raw_input("Filename for figure: ")
+    else:
+        save = False
+        fig_name = None
+
+    bw_diff_plot(file_path, var_name, colour_bounds, save, fig_name)
+
+

bwsalt_plot.py

@@ -58,7 +58,7 @@ def bwsalt_plot (file_path, save=False, fig_name=None):
     contourf(x, y, bwsalt, lev, cmap='jet', extend='both')
     cbar = colorbar(ticks=arange(min_scale, max_scale+0.2, 0.2))
     cbar.ax.tick_params(labelsize=20)
-    title(r'Bottom water salintiy (psu), annual average', fontsize=30)
+    title(r'Bottom water salinity (psu), annual average', fontsize=30)
     axis('off')
 
     # Finished

cice_grid.py

@@ -15,15 +15,38 @@ def cice_grid (roms_grid_name, cice_grid_name, cice_kmt_name):
     cice_kmt = Dataset(cice_kmt_name, 'w')
 
     # Read variables
-    ulon = roms.variables['lon_rho'][:,:]
-    ulat = roms.variables['lat_rho'][:,:]
-    angle = roms.variables['angle'][:]
+    ulon_tmp = roms.variables['lon_psi'][:,:]
+    ulat_tmp = roms.variables['lat_psi'][:,:]
+    # Convert angle to degrees
+    angle = roms.variables['angle'][:]*180/pi
     # Mask out ice shelf cavities for sea ice
     kmt = roms.variables['mask_rho'][:,:] - roms.variables['mask_zice'][:,:]
 
     i = arange(angle.shape[1])
     j = arange(angle.shape[0])
 
+    # Get one more index in each dimension for lat and lon; linearly extend
+    ulon = zeros([size(j), size(i)])
+    ulon[:-1,:-1] = ulon_tmp
+    # Longitude is a bit tricky because of the mod 360
+    for jj in range(size(j)):
+        ulon_1back = ulon[jj,-2]
+        ulon_2back = ulon[jj,-3]
+        if ulon_2back - ulon_1back > 300:
+            ulon_2back -= 360
+        ulon[jj,-1] = 2*ulon_1back - ulon_2back
+    for ii in range(size(i)):
+        ulon_1back = ulon[-2,ii]
+        ulon_2back = ulon[-3,ii]
+        if ulon_2back - ulon_1back < -300:
+            ulon_2back += 360
+        ulon[-1,ii] = 2*ulon_1back - ulon_2back        
+    ulat = zeros([size(j), size(i)])
+    ulat[:-1,:-1] = ulat_tmp
+    # Latitude is easy
+    ulat[-1,:] = 2*ulat[-2,:] - ulat[-3,:]
+    ulat[:,-1] = 2*ulat[:,-2] - ulat[:,-3]
+
     # Write variables
     cice_grid.createDimension('i', size(i))
     cice_grid.createDimension('j', size(j))

cice_strx.py

@@ -0,0 +1,74 @@
+from netCDF4 import Dataset
+from numpy import *
+from matplotlib.pyplot import *
+from rotate_vector_cice import *
+
+def cice_strx (file_path, tstep, colour_bounds=None, save=False, fig_name=None):
+
+    deg2rad = pi/180
+
+    id = Dataset(file_path, 'r')
+    lon = id.variables['ULON'][:-15,:]
+    lat = id.variables['ULAT'][:-15,:]
+    angle = id.variables['ANGLET'][:-15,:]
+    strx = id.variables['strairx'][tstep-1,:-15,:] + id.variables['strocnx'][tstep-1,:-15,:] + id.variables['strtltx'][tstep-1,:-15,:] + id.variables['strcorx'][tstep-1,:-15,:] + id.variables['strintx'][tstep-1,:-15,:]
+    stry = id.variables['strairy'][tstep-1,:-15,:] + id.variables['strocny'][tstep-1,:-15,:] + id.variables['strtlty'][tstep-1,:-15,:] + id.variables['strcory'][tstep-1,:-15,:] + id.variables['strinty'][tstep-1,:-15,:]
+    id.close()
+
+    strx_lonlat, stry_lonlat = rotate_vector_cice(strx, stry, angle)
+
+    x = -(lat+90)*cos(lon*deg2rad+pi/2)
+    y = (lat+90)*sin(lon*deg2rad+pi/2)
+
+    if colour_bounds is not None:
+        # User has set bounds on colour scale
+        lev = linspace(colour_bounds[0], colour_bounds[1], num=40)
+        if colour_bounds[0] == -colour_bounds[1]:
+            # Bounds are centered on zero, so choose a blue-to-red colourmap
+            # centered on yellow
+            colour_map = 'RdYlBu_r'
+        else:
+            colour_map = 'jet'
+    else:
+        max_val = amax(abs(stry_lonlat))
+        lev = linspace(-max_val, max_val, num=40)
+        colour_map = 'RdYlBu_r'
+
+    fig = figure(figsize=(16,12))
+    fig.add_subplot(1,1,1, aspect='equal')
+    contourf(x, y, strx_lonlat, lev, cmap=colour_map, extend='both')
+    cbar = colorbar()
+    cbar.ax.tick_params(labelsize=20)
+    title('Sum of eastward stresses N/m^2', fontsize=30)
+    axis('off')
+
+    if save:
+        fig.savefig(fig_name)
+    else:
+        fig.show()
+
+
+if __name__ == "__main__":
+
+    file_path = raw_input("Path to CICE history file: ")
+    tstep = int(raw_input("Timestep number (starting at 1): "))
+
+    # Get colour bounds if necessary
+    colour_bounds = None
+    get_bounds = raw_input("Set bounds on colour scale (y/n)? ")
+    if get_bounds == 'y':
+        lower_bound = float(raw_input("Lower bound: "))
+        upper_bound = float(raw_input("Upper bound: "))
+        colour_bounds = [lower_bound, upper_bound]
+
+    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
+
+    cice_strx(file_path, tstep, colour_bounds, save, fig_name)
+
+

cice_stry.py

@@ -0,0 +1,74 @@
+from netCDF4 import Dataset
+from numpy import *
+from matplotlib.pyplot import *
+from rotate_vector_cice import *
+
+def cice_stry (file_path, tstep, colour_bounds=None, save=False, fig_name=None):
+
+    deg2rad = pi/180
+
+    id = Dataset(file_path, 'r')
+    lon = id.variables['ULON'][:-15,:]
+    lat = id.variables['ULAT'][:-15,:]
+    angle = id.variables['ANGLET'][:-15,:]
+    strx = id.variables['strairx'][tstep-1,:-15,:] + id.variables['strocnx'][tstep-1,:-15,:] + id.variables['strtltx'][tstep-1,:-15,:] + id.variables['strcorx'][tstep-1,:-15,:] + id.variables['strintx'][tstep-1,:-15,:]
+    stry = id.variables['strairy'][tstep-1,:-15,:] + id.variables['strocny'][tstep-1,:-15,:] + id.variables['strtlty'][tstep-1,:-15,:] + id.variables['strcory'][tstep-1,:-15,:] + id.variables['strinty'][tstep-1,:-15,:]
+    id.close()
+
+    strx_lonlat, stry_lonlat = rotate_vector_cice(strx, stry, angle)
+
+    x = -(lat+90)*cos(lon*deg2rad+pi/2)
+    y = (lat+90)*sin(lon*deg2rad+pi/2)
+
+    if colour_bounds is not None:
+        # User has set bounds on colour scale
+        lev = linspace(colour_bounds[0], colour_bounds[1], num=40)
+        if colour_bounds[0] == -colour_bounds[1]:
+            # Bounds are centered on zero, so choose a blue-to-red colourmap
+            # centered on yellow
+            colour_map = 'RdYlBu_r'
+        else:
+            colour_map = 'jet'
+    else:
+        max_val = amax(abs(stry_lonlat))
+        lev = linspace(-max_val, max_val, num=40)
+        colour_map = 'RdYlBu_r'
+
+    fig = figure(figsize=(16,12))
+    fig.add_subplot(1,1,1, aspect='equal')
+    contourf(x, y, stry_lonlat, lev, cmap=colour_map, extend='both')
+    cbar = colorbar()
+    cbar.ax.tick_params(labelsize=20)
+    title('Sum of northward stresses N/m^2', fontsize=30)
+    axis('off')
+
+    if save:
+        fig.savefig(fig_name)
+    else:
+        fig.show()
+
+
+if __name__ == "__main__":
+
+    file_path = raw_input("Path to CICE history file: ")
+    tstep = int(raw_input("Timestep number (starting at 1): "))
+
+    # Get colour bounds if necessary
+    colour_bounds = None
+    get_bounds = raw_input("Set bounds on colour scale (y/n)? ")
+    if get_bounds == 'y':
+        lower_bound = float(raw_input("Lower bound: "))
+        upper_bound = float(raw_input("Upper bound: "))
+        colour_bounds = [lower_bound, upper_bound]
+
+    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
+
+    cice_stry(file_path, tstep, colour_bounds, save, fig_name)
+
+