From 2716b38a48b1707f6e65b8a2cd54821436a8ed80 Mon Sep 17 00:00:00 2001 From: Kaitlin Alexander Date: Thu, 10 Dec 2015 17:32:10 +1100 Subject: [PATCH] Added scripts to calculate ocean heat content, salt content, and TKE --- calc_z.pyc | Bin 984 -> 0 bytes file_guide.txt | 28 ++++---- plot_ohc.py | 18 +++--- plot_tke.py | 162 ++++++++++++++++++++++++++++++++++++++++++++++ plot_totalsalt.py | 18 +++--- 5 files changed, 197 insertions(+), 29 deletions(-) delete mode 100644 calc_z.pyc create mode 100644 plot_tke.py diff --git a/calc_z.pyc b/calc_z.pyc deleted file mode 100644 index bf018bf428b12cff8784f747e1b65a4c0b672e03..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 984 zcmb7B&2G~`5T13?Bu(?z#ID1ILy?fxLqiT!94ZP5=Nu?M5)z6W$HY;cIAk}4M!ToN zbMR_B0um1Z->egIX0-*3M06n@vMKl(47?NloGv}IP{I^gP=)) zu7e_Cvp+b2fTmcSMHR>syfWSZPu{|ha#OMLkfkn@XR4$SaM*zK?;Qs1nj_MCo4R zREhmCUwwQ_#${9uY?mIm8lhUQ`s^-|?-|c#c;xPT>kcRJLWsVD!A?P&%g-(Gp5wr2 zao_9U0RMkDg~0iBHRA{*~;bIzNE}|%tiB;rmk)Sn#W4U!~Vhl&L ziC3Alnk@@>;dcLd;bs=ailNhNFAu7~QMu98@VBWeA)KE}?7lE|VAEp(VJlhF1Z zBiqI&((*~;Tbq}&8p!p^rVYlKvcr6dXVXO4l!5KBSkFh`ve&~ry0JMo7JKs~00Xi@ zK8DtlNw!kG`Ngx|O&C7C==WCHT#i(hrLq@=X*AlLFYj!_-0Gl6o-#?o ziBoYVPMnt06iwlY4(z@&40v2SCb%u{UO$An@BiN 1200)*(lon_u < 100)) + lon_u[index1_u] = lon_u[index1_u] + 360 + index2_u = nonzero((i_u < 200)*(lon_u > 300)) + lon_u[index2_u] = lon_u[index2_u] - 360 + # Repeat for v grid + i_v = tile(arange(1, num_lon_v+1), (num_lat_v, 1)) + index1_v = nonzero((i_v > 1200)*(lon_v < 100)) + lon_v[index1_v] = lon_v[index1_v] + 360 + index2_v = nonzero((i_v < 200)*(lon_v > 300)) + lon_v[index2_v] = lon_v[index2_v] - 360 + + # Interpolate to get longitude at the edges of each cell + w_bdry_u = (lon_u[:,0] + lon_u[:,-1] - 360)/2 + middle_lon_u = (lon_u[:,0:-1] + lon_u[:,1:])/2 + e_bdry_u = (lon_u[:,0] + 360 + lon_u[:,-1])/2 + lon_u_edges = ma.concatenate((w_bdry_u[:,None], middle_lon_u, e_bdry_u[:,None]), axis=1) + # Subtract to get the change in longitude over each cell + dlon_u = abs(lon_u_edges[:,1:] - lon_u_edges[:,0:-1]) + # Repeat for v grid + w_bdry_v = (lon_v[:,0] + lon_v[:,-1] - 360)/2 + middle_lon_v = (lon_v[:,0:-1] + lon_v[:,1:])/2 + e_bdry_v = (lon_v[:,0] + 360 + lon_v[:,-1])/2 + lon_v_edges = ma.concatenate((w_bdry_v[:,None], middle_lon_v, e_bdry_v[:,None]), axis=1) + dlon_v = abs(lon_v_edges[:,1:] - lon_v_edges[:,0:-1]) + + # Similarly for latitude + s_bdry_u = lat_u[0,:] + middle_lat_u = (lat_u[0:-1,:] + lat_u[1:,:])/2 + n_bdry_u = lat_u[-1,:]*0 - 50 + lat_u_edges = ma.concatenate((s_bdry_u[None,:], middle_lat_u, n_bdry_u[None,:])) + dlat_u = lat_u_edges[1:,:] - lat_u_edges[0:-1,:] + # Repeat for v grid + s_bdry_v = lat_v[0,:] + middle_lat_v = (lat_v[0:-1,:] + lat_v[1:,:])/2 + n_bdry_v = lat_v[-1,:]*0 - 50 + lat_v_edges = ma.concatenate((s_bdry_v[None,:], middle_lat_v, n_bdry_v[None,:])) + dlat_v = lat_v_edges[1:,:] - lat_v_edges[0:-1,:] + + # Convert from spherical to Cartesian coordinates + # dy = r*dlat where dlat is converted to radians + dy_u_2d = r*dlat_u*pi/180.0 + dy_v_2d = r*dlat_v*pi/180.0 + # Copy into a 3D array, same at each depth level + dy_u = tile(dy_u_2d, (N,1,1)) + dy_v = tile(dy_v_2d, (N,1,1)) + # dx = r*cos(lat)*dlon where lat and dlon are converted to radians + dx_u_2d = r*cos(pi*lat_u/180.0)*dlon_u*pi/180.0 + dx_v_2d = r*cos(pi*lat_v/180.0)*dlon_v*pi/180.0 + dx_u = tile(dx_u_2d, (N,1,1)) + dx_v = tile(dx_v_2d, (N,1,1)) + + # Get a 3D array of z-coordinates on u and v grids + # sc_r and Cs_r are unused in this script + z_u, sc_r, Cs_r = calc_z(h_u, zice_u, lon_u, lat_u, theta_s, theta_b, hc, N) + z_v, sc_r, Cs_r = calc_z(h_v, zice_v, lon_v, lat_v, theta_s, theta_b, hc, N) + # We have z at the midpoint of each cell, now find it on the top and + # bottom edges of each cell + z_u_edges = zeros((size(z_u,0)+1, size(z_u,1), size(z_u,2))) + z_u_edges[1:-1,:,:] = 0.5*(z_u[0:-1,:,:] + z_u[1:,:,:]) + # At surface, z = 0; at bottom, set z to be the same as the midpoint of + # the deepest cell + z_u_edges[0,:,:] = z_u[0,:,:] + # Now find dz + dz_u = z_u_edges[1:,:,:] - z_u_edges[0:-1,:,:] + # Repeat on the v grid + z_v_edges = zeros((size(z_v,0)+1, size(z_v,1), size(z_v,2))) + z_v_edges[1:-1,:,:] = 0.5*(z_v[0:-1,:,:] + z_v[1:,:,:]) + z_v_edges[0,:,:] = z_v[0,:,:] + dz_v = z_v_edges[1:,:,:] - z_v_edges[0:-1,:,:] + # Calculate dV for each grid + dV_u = dx_u*dy_u*dz_u + dV_v = dx_v*dy_v*dz_v + + # Read time data and convert from seconds to years + id = Dataset(file_path, 'r') + time = id.variables['ocean_time'][:]/(365*24*60*60) + # Read reference density + rho0 = id.variables['rho0'][:] + + avgke = zeros(size(time)) + for t in range(size(time)): + # Read u, v, and rho at this timestep + u = id.variables['u'][t,:,:,:] + v = id.variables['v'][t,:,:,:] + rho = id.variables['rho'][t,:,:,:] + rho0 + # Interpolate rho onto u and v grids + rho_u = 0.5*(rho[:,:,0:-1] + rho[:,:,1:]) + rho_v = 0.5*(rho[:,0:-1,:] + rho[:,1:,:]) + # Integrate 0.5*rho*vel^2 over volume on each grid, add for TKE + avgke[t] = sum(0.5*rho_u*u**2*dV_u) + sum(0.5*rho_v*v**2*dV_v) + id.close() + + # Plot results + clf() + plot(time, avgke) + xlabel('Years') + ylabel('Southern Ocean Total Kinetic Energy (J)') + show() + + +# Command-line interface +if __name__ == "__main__": + + file_path = raw_input("Path to ocean averages file: ") + grid_path = raw_input("Path to grid file: ") + plot_tke(file_path, grid_path) + + + + + + diff --git a/plot_totalsalt.py b/plot_totalsalt.py index c54533c..36fe6fc 100644 --- a/plot_totalsalt.py +++ b/plot_totalsalt.py @@ -3,7 +3,7 @@ from matplotlib.pyplot import * from calc_z import * -# Calculate the total ocean salt content at each timestep of the given ocean +# Calculate and plot the total salt content at each timestep of the given ocean # averages file. def plot_totalsalt (file_path, grid_path): @@ -42,19 +42,19 @@ def plot_totalsalt (file_path, grid_path): lon[index2] = lon[index2] - 360 # Interpolate to get longitude at the edges of each cell - w_bdry = (lon[:,0] + lon[:,num_lon-1] - 360)/2 - middle_lon = (lon[:,0:num_lon-1] + lon[:,1:num_lon])/2 - e_bdry = (lon[:,0] + 360 + lon[:,num_lon-1])/2 + w_bdry = (lon[:,0] + lon[:,-1] - 360)/2 + middle_lon = (lon[:,0:-1] + lon[:,1:])/2 + e_bdry = (lon[:,0] + 360 + lon[:,-1])/2 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:num_lon+1] - lon_edges[:,0:num_lon]) + dlon = abs(lon_edges[:,1:] - lon_edges[:,0:-1]) # Similarly for latitude s_bdry = lat[0,:] - middle_lat = (lat[0:num_lat-1,:] + lat[1:num_lat,:])/2 - n_bdry = lat[num_lat-1,:]*0 - 50 + middle_lat = (lat[0:-1,:] + lat[1:,:])/2 + n_bdry = lat[-1,:]*0 - 50 lat_edges = ma.concatenate((s_bdry[None,:], middle_lat, n_bdry[None,:])) - dlat = lat_edges[1:num_lat+1,:] - lat_edges[0:num_lat,:] + dlat = lat_edges[1:,:] - lat_edges[0:-1,:] # Convert from spherical to Cartesian coordinates # dy = r*dlat where dlat is converted to radians @@ -97,7 +97,7 @@ def plot_totalsalt (file_path, grid_path): clf() plot(time, totalsalt) xlabel('Years') - ylabel('Ocean Salt Content (kg)') + ylabel('Southern Ocean Salt Content (kg)') show()