Added overturning streamfunction plot; also a few bug fixes

circumpolar_cice_plot.py

@@ -45,7 +45,7 @@ def circumpolar_cice_plot (file_path, var_name, tstep, save=False, fig_name=None
     y = (lat+90)*sin(lon*deg2rad+pi/2)
 
     # Center levels on 0 for certain variables, with a blue-to-red colourmap
-    if var_name in ['uvel', 'vvel', 'uatm', 'vatm', 'uocn', 'vocn']:
+    if var_name in ['uvel', 'vvel', 'uatm', 'vatm', 'uocn', 'vocn', 'frzmlt']:
         max_val = amax(abs(data))
         lev = linspace(-max_val, max_val, num=40)
         colour_map = 'RdYlBu_r'
@@ -66,7 +66,7 @@ def circumpolar_cice_plot (file_path, var_name, tstep, save=False, fig_name=None
         savefig(fig_name)
         close()
     else:
-        show()
+        show(block=False)
 
 
 # Command-line interface

circumpolar_plot.py

@@ -156,7 +156,7 @@ def circumpolar_plot (file_path, var_name, tstep, depth_key, depth, depth_bounds
         savefig(fig_name)
         close()
     else:
-        show()
+        show(block=False)
 
 
 # Linearly interpolate data to the specified depth at each horizontal point.

file_guide.txt

@@ -463,6 +463,15 @@ cmip5_all_plots.py: Call cmip5_plot.py for all models (including the multi-model
 		            "run cmip5_all_plots.py".
 
 
+overturning_plot.py: Calculate the meridional overturning streamfunction at the
+                     last timestep of the given ROMS history/averages file and
+		     make a contour plot in z-space.
+		     To run: Open python or ipython and type
+		             "run overturning_plot.py". You will be prompted
+			     for the path to the ROMS history/averages file
+			     and the filename with which to save the plot.
+
+
 ***OBSOLETE***
 
 roms_grid_rtopo2.py: Interpolate RTopo-2 data to an existing ROMS grid, and do

ice2ocn_fwflux.py

@@ -35,7 +35,8 @@ def ice2ocn_fwflux (file_path, tstep, fig_name):
     title('CICE-to-ROMS net freshwater flux (cm/day)', fontsize=30)
     axis('off')
 
-    savefig(fig_name)
+    #savefig(fig_name)
+    show()
 
 
 # Command-line interface

ini_sst_circumpolar.py

@@ -49,7 +49,8 @@ def ini_sst_circumpolar (grid_path, file_path, fig_name):
     title(r'Initial SST ($^{\circ}$C)', fontsize=30)
     axis('off')
 
-    savefig(fig_name)
+    #savefig(fig_name)
+    show()
 
 
 # Command-line interface

overturning_plot.py

@@ -0,0 +1,107 @@
+from netCDF4 import Dataset
+from numpy import *
+from matplotlib.pyplot import *
+from calc_z import *
+
+# Calculate the meridional overturning streamfunction at the last timestep
+# of the given ROMS history/averages file and make a contour plot in z-space.
+# Input:
+# file_path = path to ROMS history/averages file
+# fig_name = filename for figure 
+def overturning_plot (file_path, fig_name):
+
+    # Grid parameters
+    theta_s = 0.9
+    theta_b = 4.0
+    hc = 40
+    N = 31
+    # Radius of the Earth in metres
+    r = 6.371e6
+    # Degrees to radians conversion factor
+    deg2rad = pi/180.0
+
+    # Read v and grid variables
+    id = Dataset(file_path, 'r')
+    v = id.variables['v'][-1,:,:,:]  # -1 means last timestep
+    h = id.variables['h'][:,:]
+    zice = id.variables['zice'][:,:]
+    zeta = id.variables['zeta'][-1,:,:]
+    lon = id.variables['lon_v'][:,:]
+    lat = id.variables['lat_v'][:,:]
+    id.close()
+
+    num_lat = size(lon,0)
+    num_lon = size(lon,1)
+
+    # Add or subtract 360 from longitude values which wrap around
+    # so that longitude increases monotonically from west to east
+    i = tile(arange(1,num_lon+1), (num_lat,1))
+    index1 = nonzero((i > 1200)*(lon < 100))
+    lon[index1] = lon[index1] + 360
+    index2 = nonzero((i < 200)*(lon > 300))
+    lon[index2] = lon[index2] - 360
+
+    # Interpolate to get longitude at the edges of each cell
+    w_bdry = 0.5*(lon[:,0] + lon[:,-1] - 360)
+    middle_lon = 0.5*(lon[:,0:-1] + lon[:,1:])
+    e_bdry = 0.5*(lon[:,0] + 360 + lon[:,-1])
+    lon_edges = ma.concatenate((w_bdry[:,None], middle_lon, e_bdry[:,None]), axis=1)
+    # Subtract to get the change in longitude over each cell
+    dlon = abs(lon_edges[:,1:] - lon_edges[:,0:-1])
+    # dx = r*cos(lat)*dlon where lat and dlon are converted to radians
+    dx_2d = r*cos(lat*deg2rad)*dlon*deg2rad
+    # Copy into a 3D array of dimension depth x latitude x longitude
+    dx = tile(dx_2d, (N,1,1))
+
+    # Interpolate h, zice, and zeta to the v-grid
+    h_v = 0.5*(h[0:-1,:] + h[1:,:])
+    zice_v = 0.5*(zice[0:-1,:] + zice[1:,:])
+    zeta_v = 0.5*(zeta[0:-1,:] + zeta[1:,:])
+
+    # Get a 3D grid of z-coordinates; sc_r and Cs_r are unused in this script
+    z, sc_r, Cs_r = calc_z(h_v, zice_v, theta_s, theta_b, hc, N, zeta_v)
+    # We have z at the midpoint of each cell, now find it on the top and
+    # bottom edges of each cell
+    z_edges = zeros((size(z,0)+1, size(z,1), size(z,2)))
+    z_edges[1:-1,:,:] = 0.5*(z[0:-1,:,:] + z[1:,:,:])
+    # At the surface, z = zice; at the bottom, extrapolate
+    z_edges[-1,:,:] = zice_v[:,:]
+    z_edges[0,:,:] = 2*z[0,:,:] - z_edges[1,:,:]
+    # Now find dz
+    dz = z_edges[1:,:,:] - z_edges[0:-1,:,:]
+
+    # Calculate transport in each cell
+    transport = v*dx*dz
+    # Definite integral over longitude
+    transport = sum(transport, axis=2)
+    # Indefinite integral over depth; flip before and after so the integral
+    # starts at the surface, not the bottom. Also convert to Sv.
+    transport = flipud(cumsum(flipud(transport), axis=0))*1e-6
+
+    # Calculate latitude and z coordinates, averaged over longitude,
+    # for plotting
+    avg_lat = mean(lat, axis=1)
+    avg_lat = tile(avg_lat, (N,1))
+    avg_z = mean(z, axis=2)
+
+    # Centre colour scale on 0
+    max_val = amax(abs(transport))
+    lev = linspace(-max_val, max_val, num=40)
+
+    # Make the plot
+    figure(figsize=(16,8))
+    contourf(avg_lat, avg_z, transport, lev, cmap='RdBu_r')
+    colorbar()
+    xlabel('Latitude')
+    ylabel('Depth (m)')
+    title('Meridional Overturning Streamfunction (Sv)')
+
+    savefig(fig_name)
+
+
+# Command-line interface
+if __name__ == "__main__":
+
+    file_path = raw_input("Path to ROMS history/averages file: ")
+    fig_name = raw_input("File name for figure: ")
+    overturning_plot(file_path, fig_name)

romscice_ini_ecco.py

@@ -20,6 +20,13 @@
 # Main routine
 def run (grid_file, theta_file, salt_file, output_file, Tcline, theta_s, theta_b, hc, N, nbdry_ecco):
 
+    # Parameters for freezing point of seawater; taken from Ben's PhD thesis
+    a = -5.73e-2
+    b = 8.32e-2
+    c = -7.61e-8
+    rho = 1035.0
+    g = 9.81    
+
     # Read ECCO2 data and grid
     print 'Reading ECCO2 data'
     theta_fid = Dataset(theta_file, 'r')
@@ -92,7 +99,20 @@ def run (grid_file, theta_file, salt_file, output_file, Tcline, theta_s, theta_b
     print 'Interpolating temperature'
     temp = interp_ecco2roms(theta, lon_ecco, lat_ecco, depth_ecco, lon_roms_3d, lat_roms_3d, z_roms_3d, -0.5)
     print 'Interpolating salinity'
-    salt = interp_ecco2roms(salt, lon_ecco, lat_ecco, depth_ecco, lon_roms_3d, lat_roms_3d, z_roms_3d, 34.5)
+    salt = interp_ecco2roms(salt, lon_ecco, lat_ecco, depth_ecco, lon_roms_3d, lat_roms_3d, z_roms_3d, 35)
+
+    # Identify the continental shelf
+    h_3d = tile(h, (N,1,1))
+    index = nonzero((lat_roms_3d <= -60)*(h_3d <= 1200))
+    # Set salinity over the continental shelf to a constant value, and
+    # temperature to the local freezing point
+    salt[index] = 35
+    temp[index] = a*salt[index] + b + c*rho*g*abs(z_roms_3d[index])
+    # Make sure the ice shelf cavities got this too
+    mask_zice_3d = tile(mask_zice, (N,1,1))
+    index = mask_zice_3d == 1
+    salt[index] = 35
+    temp[index] = a*salt[index] + b + c*rho*g*abs(z_roms_3d[index])
 
     # Set initial velocities and sea surface height to zero
     u = zeros((N, num_lat, num_lon-1))
@@ -186,8 +206,8 @@ def run (grid_file, theta_file, salt_file, output_file, Tcline, theta_s, theta_b
     print 'Finished'
 
 
-# Given an array on the ECCO2 grid, fill in the land mask with constant values,
-# and then interpolate to the ROMS grid.
+# Given an array on the ECCO2 grid, fill land mask with constant values
+# and interpolate onto the ROMS grid.
 # Input:
 # A = array of size nxmxo containing values on the ECCO2 grid (dimension
 #     longitude x latitude x depth)
@@ -203,10 +223,8 @@ def run (grid_file, theta_file, salt_file, output_file, Tcline, theta_s, theta_b
 # z_roms_3d = array of size pxqxr containing ROMS depth values in negative
 #             z-coordinates (converted from grid file using calc_z.py, as
 #             shown below in the main script)
-# fill_value = scalar containing the value with which to fill the ECCO2 land
-#              mask: something close to the mean value of A so that the
-#              splines don't freak out (suggest -0.5 for temperature and 34.5
-#              for salinity)
+# fill_value = scalar containing the value with which to fill the ROMS land mask
+#              and ice shelf cavities (doesn't really matter, don't use NaN)
 # Output:
 # B = array of size pxqxr containing values on the ROMS grid (dimension depth x
 #     latitude x longitude)
@@ -215,24 +233,21 @@ def interp_ecco2roms (A, lon_ecco, lat_ecco, depth_ecco, lon_roms_3d, lat_roms_3
     # Calculate N based on size of ROMS grid
     N = size(lon_roms_3d, 0)
 
-    # Fill the ECCO2 land mask with a constant value close to the mean of A.
-    # This is less than ideal because it will skew the interpolation slightly
-    # along the coast, and the land masks won't exactly line up between the
-    # grids. However it allows us to use the MUCH faster RegularGridInterpolator
-    # (which requires data on a regular grid i.e. no missing values) and this
-    # initialisation file is designed for a partial spinup after all!
+    # Fill the ECCO2 land mask with a constant value close to the mean of A
+    # This is less than ideal, but we will overwrite the continental shelf
+    # values anyway later
     Afill = A
     Afill[A.mask] = fill_value
 
-    # Build a function for linear interpolation of Afill
+    # Build a function for linear interpolation of A
     interp_function = RegularGridInterpolator((lon_ecco, lat_ecco, depth_ecco), Afill)
     B = zeros(shape(lon_roms_3d))
-    # Call this function once at each grid level - 3D vectorisation uses too
+    # Interpolate each z-level individually - 3D vectorisation uses too
     # much memory!
     for k in range(N):
         print '...vertical level ', str(k+1), ' of ', str(N)
         # Pass positive values for ROMS depth
-        B[k,:,:] = interp_function((lon_roms_3d[k,:,:], lat_roms_3d[k,:,:], -z_roms_3d[k,:,:]))    
+        B[k,:,:] = interp_function((lon_roms_3d[k,:,:], lat_roms_3d[k,:,:], -z_roms_3d[k,:,:]))
 
     return B
 
@@ -242,12 +257,12 @@ def interp_ecco2roms (A, lon_ecco, lat_ecco, depth_ecco, lon_roms_3d, lat_roms_3
 
     # File paths to edit here:
     # Path to ROMS grid file
-    grid_file = '/short/m68/kaa561/roms_circumpolar/data/caisom001_OneQuartergrd.nc'
+    grid_file = '../ROMS-CICE-MCT/apps/common/grid/circ30S_quarterdegree_rp5.nc'
     # Path to ECCO2 files for temperature and salinity in January 1995
-    theta_file = '../data/ECCO2/THETA.1440x720x50.199501.nc'
-    salt_file = '../data/ECCO2/SALT.1440x720x50.199501.nc'
+    theta_file = '../ROMS-CICE-MCT/data/ECCO2/raw/THETA.1440x720x50.199501.nc'
+    salt_file = '../ROMS-CICE-MCT/data/ECCO2/raw/SALT.1440x720x50.199501.nc'
     # Path to desired output file
-    output_file = '../data/caisom001_ini_1995.nc'
+    output_file = '../ROMS-CICE-MCT/data/ecco2_ini.nc'
 
     # Grid parameters; check grid_file and *.in file to make sure these are correct
     Tcline = 40

romscice_ini_woa.py

@@ -63,7 +63,7 @@ def run (grid_file, woa_file, output_file, Tcline, theta_s, theta_b, hc, N, nbdr
     print 'Interpolating temperature'
     temp = interp_woa2roms(temp_woa, lon_woa, lat_woa, depth_woa, lon_roms_3d, lat_roms_3d, z_roms_3d, mask_rho, mask_zice, -0.5)
     print 'Interpolating salinity'
-    salt = interp_woa2roms(salt_woa, lon_woa, lat_woa, depth_woa, lon_roms_3d, lat_roms_3d, z_roms_3d, mask_rho, mask_zice, 34.5)
+    salt = interp_woa2roms(salt_woa, lon_woa, lat_woa, depth_woa, lon_roms_3d, lat_roms_3d, z_roms_3d, mask_rho, mask_zice, 35)
 
     # Set the temperature in ice shelf cavities to the local freezing point
     mask_zice_3d = tile(mask_zice, (N,1,1))

spinup_plots.py

@@ -120,8 +120,7 @@ def calc_grid (file_path):
 def get_rho (file_path, t):
 
     # Reference density
-    # New users: edit this based on the value in your .in ROMS configuration file
-    rho0 = 1025.0
+    rho0 = 1000.0
 
     id = Dataset(file_path, 'r')
     # Read density anomalies, add rho0 to get absolute density

zonal_plot.py

@@ -169,7 +169,7 @@ def zonal_plot (file_path, var_name, tstep, lon_key, lon0, lon_bounds, depth_min
         savefig(fig_name)
         close()
     else:
-        show()
+        show(block=False)
 
     # Reset lon0 or lon_bounds to (-180, 180) range in case we
     # use them again for the next plot